search

liucongcas / GLORI-tools

bioinformatic pipeline for GLORI
MIT License
17 stars 3 forks source link

Empty _tf_rs.bam file #13

Open yuanhao96 opened 1 year ago

yuanhao96 commented 1 year ago

Thank you for developing the pipeline! I encountered an error at the alignment step of runGLORI.py. It seems that one of the mapping results (_tf_rs.bam) is empty. The reference genome and fastq files were preprocessed following the README file, except that base-annotation was skipped.

Could you please give me some hints what may cause this error, or is there any walkaround to address the problem?

Please find below my log before returning the error:

**************combine,treated
python /nfs/turbo/umms-rajeshr/softwares/GLORI-tools/pipelines/mapping_reads.py -q /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/preprocessing/fastx_trimmer/8324-JX-1_CTTACCGA-GATTCCGA_S28_R2_001_trimmed.fq -p 30 -f /nfs/turbo/umms-rajeshr/reference_genome/refdata-cellranger-mm10-3.0.0/mm10_GLORI/mm10.AG_conversion.fa --FilterN 0.5 -rvs /nfs/turbo/umms-rajeshr/reference_genome/refdata-cellranger-mm10-3.0.0/mm10_GLORI/mm10.rvsCom.fa -Tf /nfs/turbo/umms-rajeshr/reference_genome/refdata-cellranger-mm10-3.0.0/mm10_GLORI/mm10.cdna2.fa.AG_conversion.fa -mulMax 1 -t STAR -m 2 -pre 8324-JX-1 -o /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1 --combine  --rvs_fac
    STAR --runThreadN 30 --genomeDir /nfs/turbo/umms-rajeshr/reference_genome/refdata-cellranger-mm10-3.0.0/mm10_GLORI/mm10.AG_conversion --limitOutSJcollapsed 5000000 --outFilterMismatchNmax 2 --outFilterScoreMinOverLread 0.5 --outFilterMatchNminOverLread 0.5 --seedSearchStartLmax 30 --outSAMattributes All --outSAMprimaryFlag AllBestScore --outMultimapperOrder Random --outSAMmultNmax 1 --outSAMtype BAM Unsorted --outFilterMultimapNmax 1 --outFileNamePrefix /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1. --readFilesIn /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1//8324-JX-1_AGchanged_2.fq --outSAMunmapped Within --outReadsUnmapped Fastx
    STAR version: 2.7.10b   compiled: 2022-11-01T09:53:26-04:00 :/home/dobin/data/STAR/STARcode/STAR.master/source
Jul 11 11:06:23 ..... started STAR run
Jul 11 11:06:23 ..... loading genome
Jul 11 11:07:34 ..... started mapping
Jul 11 11:20:54 ..... finished mapping
Jul 11 11:20:57 ..... finished successfully
[bam_sort_core] merging from 12 files and 1 in-memory blocks...
[bam_sort_core] merging from 12 files and 1 in-memory blocks...
[bam_sort_core] merging from 12 files and 1 in-memory blocks...
    STAR --runThreadN 30 --genomeDir /nfs/turbo/umms-rajeshr/reference_genome/refdata-cellranger-mm10-3.0.0/mm10_GLORI/mm10.rvsCom --limitOutSJcollapsed 5000000 --outFilterMismatchNmax 2 --outFilterScoreMinOverLread 0.5 --outFilterMatchNminOverLread 0.5 --seedSearchStartLmax 30 --outSAMattributes All --outSAMprimaryFlag AllBestScore --outMultimapperOrder Random --outSAMmultNmax 1 --outSAMtype BAM Unsorted --outFilterMultimapNmax 1 --outFileNamePrefix /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_rvs. --readFilesIn /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_un_2.fq --outSAMunmapped Within --outReadsUnmapped Fastx
    STAR version: 2.7.10b   compiled: 2022-11-01T09:53:26-04:00 :/home/dobin/data/STAR/STARcode/STAR.master/source
Jul 11 12:47:47 ..... started STAR run
Jul 11 12:47:47 ..... loading genome
Jul 11 12:49:00 ..... started mapping
Jul 11 12:59:12 ..... finished mapping
Jul 11 12:59:15 ..... finished successfully
[bam_sort_core] merging from 7 files and 1 in-memory blocks...
[bam_sort_core] merging from 7 files and 1 in-memory blocks...
[bam_sort_core] merging from 7 files and 1 in-memory blocks...
[E::hts_open_format] Failed to open file /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_tf.sam
samtools view: failed to open "/nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_tf.sam" for reading: No such file or directory
[2023-07-11 11:06:22] Warning：the changed files already exists, please make sure the input file is correct
[2023-07-11 11:06:22] map to genome...
[2023-07-11 12:47:47]map to reversegenome...
[2023-07-11 13:53:45]map to transcriptome...
STAR --runThreadN 30 --genomeDir /nfs/turbo/umms-rajeshr/reference_genome/refdata-cellranger-mm10-3.0.0/mm10_GLORI/mm10.AG_conversion --limitOutSJcollapsed 5000000  --outFilterMismatchNmax 2 --outFilterScoreMinOverLread 0.5 --outFilterMatchNminOverLread 0.5 --seedSearchStartLmax 30  --outSAMattributes All --outSAMprimaryFlag AllBestScore --outMultimapperOrder Random --outSAMmultNmax 1 --outSAMtype BAM Unsorted --outFilterMultimapNmax 1 --outFileNamePrefix /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1. --readFilesIn /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1//8324-JX-1_AGchanged_2.fq --outSAMunmapped Within --outReadsUnmapped Fastx
samtools view -F 20 -@ 30 -h /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1.Aligned.out.bam | samtools sort -n -O SAM > /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1.sam
samtools view -F 20 -@ 30 -h /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1.sam | samtools sort -n -O SAM > /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_sorted.sam
/nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_sorted.sam /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1//8324-JX-1_A.bed_sorted
************reversed_reads == mapped reads*************** 26845062 26845062
samtools view -F 20 -bS -@ 30 -h /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_r.sam | samtools sort > /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_rs.bam
STAR --runThreadN 30 --genomeDir /nfs/turbo/umms-rajeshr/reference_genome/refdata-cellranger-mm10-3.0.0/mm10_GLORI/mm10.rvsCom --limitOutSJcollapsed 5000000  --outFilterMismatchNmax 2 --outFilterScoreMinOverLread 0.5 --outFilterMatchNminOverLread 0.5 --seedSearchStartLmax 30  --outSAMattributes All --outSAMprimaryFlag AllBestScore --outMultimapperOrder Random --outSAMmultNmax 1 --outSAMtype BAM Unsorted --outFilterMultimapNmax 1 --outFileNamePrefix /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_rvs. --readFilesIn /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_un_2.fq --outSAMunmapped Within --outReadsUnmapped Fastx
samtools view -F 4 -@ 30 -h /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_rvs.Aligned.out.bam | samtools sort -n -O SAM > /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_rvs.sam
samtools view -F 4 -q 255 -@ 30 -h /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_rvs.sam | samtools sort -n -O SAM > /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_rvs_sorted.sam
/nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_rvs_sorted.sam /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1//8324-JX-1_A.bed_sorted
************reversed_reads == mapped reads*************** 16547386 16547386
samtools view -F 4 -bS -@ 30 -h /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_rvs_r.sam | samtools sort > /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_rvs_rs.bam
bowtie -k 1 -m 1 -v 2 --best --strata -p 30 -x /nfs/turbo/umms-rajeshr/reference_genome/refdata-cellranger-mm10-3.0.0/mm10_GLORI/mm10.cdna2.fa.AG_conversion.fa /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_rvs_un_2.fq -S /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_tf.sam --un /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_tf_un_2.fq 2>/nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_tf.sam.output
samtools view -F 20 -@ 30 -h /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_tf.sam | samtools sort -n -O SAM > /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_tf_sorted.sam
/nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_tf_sorted.sam /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1//8324-JX-1_A.bed_sorted
************reversed_reads == mapped reads*************** 1 1
samtools view -F 20 -bS -@ 30 -h /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_tf_r.sam | samtools sort > /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_tf_rs.bam
python /nfs/turbo/umms-rajeshr/softwares/GLORI-tools/pipelines/Transcrip2genome.py --input /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1_tf_rs.bam --output /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1.trans2Genome.bam --anno /nfs/turbo/umms-rajeshr/reference_genome/refdata-cellranger-mm10-3.0.0/mm10_GLORI/genes.gtf.tbl2 --fasta /nfs/turbo/umms-rajeshr/reference_genome/refdata-cellranger-mm10-3.0.0/mm10_GLORI/mm10.AG_conversion.fa --sort --index
[2023-07-11 13:53:49]Loading annotations...
[2023-07-11 13:53:52]Buidling genome header...
Traceback (most recent call last):
  File "/nfs/turbo/umms-rajeshr/softwares/GLORI-tools/pipelines/Transcrip2genome.py", line 355, in <module>
    with pysam.AlignmentFile(options.input, in_mode) as INPUT, pysam.AlignmentFile(options.output, out_mode,header=header) as OUTPUT,\
  File "pysam/libcalignmentfile.pyx", line 741, in pysam.libcalignmentfile.AlignmentFile.__cinit__
  File "pysam/libcalignmentfile.pyx", line 990, in pysam.libcalignmentfile.AlignmentFile._open
ValueError: file has no sequences defined (mode='rb') - is it SAM/BAM format? Consider opening with check_sq=False
[2023-07-11 13:58:01]Sorting bam...
[2023-07-11 13:58:44]Indexing bam...
[E::hts_open_format] Failed to open file /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1.trans2Genome.sorted.bam
Traceback (most recent call last):
  File "/nfs/turbo/umms-rajeshr/softwares/GLORI-tools/pipelines/concat_bam.py", line 72, in <module>
    hid,new_header,hid_dict,lift_over = read_headers(fn,hid,new_header,hid_dict,lift_over)
  File "/nfs/turbo/umms-rajeshr/softwares/GLORI-tools/pipelines/concat_bam.py", line 35, in read_headers
    with pysam.AlignmentFile(fn, 'rb') as INPUT:
  File "pysam/libcalignmentfile.pyx", line 741, in pysam.libcalignmentfile.AlignmentFile.__cinit__
  File "pysam/libcalignmentfile.pyx", line 940, in pysam.libcalignmentfile.AlignmentFile._open
FileNotFoundError: [Errno 2] could not open alignment file `/nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1.trans2Genome.sorted.bam`: No such file or directory
rm -f /nfs/turbo/umms-rajeshr/proj/jing_m6a_GRBD-Mettl3/processed/glori/8324-JX-1/8324-JX-1.trans2Genome.sorted.bam*
mmyoung commented 1 year ago

I also ran into this problem, debugging right now... Have you solved this already?

mmyoung commented 1 year ago

I've found the reason why I didn't get the _tf_rs.bam file, it's because the parameters -m --best --strata are not proper for bowtie2 in mapping_files function in pipelines/mapping_reads.py script. You may need to modify a bit of those parameters.