Rsubread: «ОШИБКА: невозможно завершить файл SAM!» - PullRequest
Новая
       20

Rsubread: «ОШИБКА: невозможно завершить файл SAM!»

0 голосов
/ 07 ноября 2018

Я хочу получить количество считываний для данных RNA-Seq, используя пакет Rsubread. Я получил свои данные RNA-Seq от sra и все другие файлы (гены, геном, ...) из сборки NCBI. Я попробовал следующий код на 4 разных организмах, и для 2 он работал отлично, но для двух других я получил следующую ошибку:

ОШИБКА: невозможно завершить файл SAM! Пожалуйста, проверьте дисковое пространство в выходном каталоге. Выходной файл не был создан.
Fehler в .load.delete.summary (output_file [i]): ОШИБКА: файл сводки Bacillus_amyloliquefaciens_samfile.sam.summary не был создан! Программа остановлена ​​неправильно!

На диске достаточно свободного места. Вот мой код: 

rna_counts <- function(path, organism, genome_gff, genome_fna, feature_table, reads_1, reads_2=NULL){
  require(Rsubread)
  require(rtracklayer)
  require(biomaRt)
  RNAseqDATADIR <- path
  setwd(RNAseqDATADIR)
  dir(RNAseqDATADIR)

  REF_GENOME <- paste0("ref_data/",genome_fna)
  RSUBREAD_INDEX_PATH <- "ref_data"
  RSUBREAD_INDEX_BASE <- organism
  buildindex(basename=file.path(RSUBREAD_INDEX_PATH,RSUBREAD_INDEX_BASE), reference=REF_GENOME)

  inputfilefwd <- file.path(RNAseqDATADIR, paste0("raw_data/", reads_1))
  if (!is.na(reads_2)) inputfilervs <- file.path(RNAseqDATADIR, paste0("/raw_data/", reads_2))
  else inputfilervs <- NA

  align(index=file.path(RSUBREAD_INDEX_PATH,RSUBREAD_INDEX_BASE), readfile1=inputfilefwd, readfile2=inputfilervs, output_file=paste0(organism,"_samfile.sam"), output_format="SAM",nthreads = 4)

  outputsamfile <- paste0(path,"/", organism, "_samfile.sam")
  propmap <- propmapped(outputsamfile)

  feature_table <- read.table(paste0(path,"/ref_data/",feature_table),header=TRUE,sep = c("\t","\n"))
  GENOME_GFF <- readGFF(paste0(path,"/ref_data/",genome_gff))

  saf_file = gff_to_saf(GENOME_GFF)       # own function for creating a saf file
  if is.na(reads_2)) paired = FALSE else paired = TRUE
  feature_counts <-featureCounts(outputsamfile, annot.ext=saf_file, isPairedEnd=paired)
  return(feature_counts)
}

Введите: 

genome_gff  = "GCA_000221645.1_ASM22164v1_genomic.gff"  
genome_fna = "GCA_000221645.1_ASM22164v1_genomic.fna"  
feature_table = GCA_000221645.1_ASM22164v1_feature_table.txt"  
reads_1 = "SRR1916792_1.fastq"  
reads_2 = "SRR1916792_2.fastq"

Во время расчетов были также распечатки, которые выглядят немного глючно: 

|| 311% completed, 17 mins elapsed, rate=7,8k fragments per second ||
|| 363% completed, 17 mins elapsed, rate=8,8k fragments per second ||
||365% completed, 18 mins elapsed, rate=8,7k fragments per second ||

Спасибо за вашу помощь!

...