Open ljudevitluka opened 4 years ago
alexdobin
commented
4 years ago Hi @ljudevitluka
the negative alignment coordinates are definitely wrong. Please try to map this suspicious read with --outSAMtype SAM option and post the result.
Cheers Alex
ljudevitluka
commented
4 years ago @alexdobin Thank you for your answer. Here are the results of mapping with --outSAMtype SAM flag. Mapping results are the some: Number of input reads | 47504246 Average input read length | 294 UNIQUE READS: Uniquely mapped reads number | 35178728 Uniquely mapped reads % | 74.05% Average mapped length | 283.67 Number of splices: Total | 17790250 Number of splices: Annotated (sjdb) | 0 Number of splices: GT/AG | 17461051 Number of splices: GC/AG | 243143 Number of splices: AT/AC | 2134 Number of splices: Non-canonical | 83922 Mismatch rate per base, % | 0.48% Deletion rate per base | 0.03% Deletion average length | 3.10 Insertion rate per base | 0.02% Insertion average length | 2.26 MULTI-MAPPING READS: Number of reads mapped to multiple loci | 7402594 % of reads mapped to multiple loci | 15.58% Number of reads mapped to too many loci | 107642 % of reads mapped to too many loci | 0.23% Again multiple reads have negative coordinates. A00709:52:H2KFFDSXY:4:1101:4291:1000 103 -1000802508 255 1S124M413N24M -1000802401 661 ACAGAGTGCCGCTGAAGAAGGTGGAGCGGAAGAGGACCCTCCAGGACCTCCGTCGCTCCCACGCCATCCTGGCCCATCGCTACGGCGCCAAGGACGACATCATCGTCATCGAGGACTACCAAGACGCTCAGTACTATGGACCCATTAGC F,F,FFFF:FFFF:FF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,FF:FFFFFFFFFFFFFFFFFFFF:FFFFFF,F:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF NH:i:1 HI:i:1 AS:i:285 nM:i:1 XS:A:+ A00709:52:H2KFFDSXY:4:1101:4291:1000 151 -1000802401 255 17M413N124M9S -1000802508 -661 TCGAGGACTACCAAGACGCTCAGTACTATGGACCCATTAGCATCGGGACGCCCGGGCAGGGGTTCGATATCATCTTCGACACTGGATCGTCCAACCTTTGGGTGCCGTCTGAACAGTGCAGCATCATCGACCTGGCCTGTCAGACGCACG FFFFFFFFF:FFF:FFFFFFF:FFFFFFFFFFFFFF:F,FF,FFFFFFFFFFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFF:FF:F:F,FFFFFFFFFF,FFFFFFFFF:FFFFFFFFFFFFFFFF NH:i:1 HI:i:1 AS:i:285 nM:i:1 XS:A:+ A00709:52:H2KFFDSXY:4:1101:5683:1000 99 SEQ13 178295324 255 47M2I99M = 178295511 328 AAATACCATATACTGTACTTAAATGTATGAAAAAGTTGAATATTATTTAAAAAAAAATAGAGTAGAATTCTTTTTTGAATGTGCTGTAACTGTTTGTTTCATTAATATACTGAAAGATACATTTCAGTGAAGTGTAGTGATAAGAGAG FFFFFF:FFFFFFF:FFFF,FFFFFFFFFFFF:FFFFFFFFF:FFFFFFFFFFFFFFFFFFFF,FFFFF:FFFFFFFFFFFFFFFFFFFFF:FF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF NH:i:1 HI:i:1 AS:i:269 nM:i:5 A00709:52:H2KFFDSXY:4:1101:5683:1000 147 SEQ13 178295511 255 141M9S = 178295324 -328 TTCACTCATTTTAGGCTGATTTCCTATTTATAGAAACCTATTGTATTTTAAAACTACTGTAGTCTACTAGAACATTATCTCAGAAGAGAAGAGATGTCCCATGTGACAGCTCACTTAAGGTTAAAACCATATTTATTTAGTTTAGTGATT FFFFF:FFFF:FF:F:FFF,FFF:FFFFFFFFFFFF,,FFFF:FFFFF:FFF:FFFFFFFFFF,FF,FFFFFFFFFFFFFFFFFFFFFF,FFFF:,FFFFFFFF:FFFFFFFFF:FFFF:FFFFFFFFFFF:FFFFFFFFFFFFFFFFFF NH:i:1 HI:i:1 AS:i:269 nM:i:5 A00709:52:H2KFFDSXY:4:1101:6316:1000 99 SEQ13 1093683801 255 1S148M = 1093683858 207 TTGTATTTCATCAACAACATCCAGCCAGTTGGAGAGACGCTCAGCGTCAAAGCGAGCAGTAAGTTGGTGGTGAGCCCAGAAGAAGAGCTCTCCCTTGCGGTCCAGGTGGTAGCCGTAGGAGTCCTTCCACCAGAAGGGGAAGTCCATAT FFFFFFFFFFFFFFFFFFFFFFFFFF:F,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,FFFFFFFF NH:i:1 HI:i:1 AS:i:290 nM:i:3 A00709:52:H2KFFDSXY:4:1101:6316:1000 147 SEQ13 1093683858 255 150M = 1093683801 -207 GTAAGTTGGTGGTGAGCCCAGAAGAAGAGCTCTCCCTTGCGGTCCAGGTGGTAGCCGTAGGAGTCCTTCCACCAGAAGGGGAAGTCCATATGCCAAGTAACGTGGTGAATGTTCATGCCAATATCCTCGCCAAAGTAGGCCACTCGCTGT FFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFF,FFF:,FFFF:FFFFFFFF:FFFFFFFFFFFFFF,FFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFF NH:i:1 HI:i:1 AS:i:290 nM:i:3 A00709:52:H2KFFDSXY:4:1101:16532:1000 167 -997713677 3 150M -997713494 333 AATTACTGAAAGATAGAAACCAACCTGGCTTACACCGGTCTGAACTCAAATCATGTAAAATTTTAAAGGTCGAACAGACCTTCTATTATAGTCTCTACTCTATAAAGAACTTTTAATTCAACATCGAGGTCGCAAACTTTCTTGTCGATA :FFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFF:FFFF:FFFFFFFF,FFFFFF:FFFFFFFFFF,FFFFFFF:FFFF:,FFFFFFFFFFFFFFFFFFFFF:FFFFFF:F:FFFFF:FF::,F,:F:FFFFFFFFFFFF NH:i:2 HI:i:1 AS:i:298 nM:i:0 A00709:52:H2KFFDSXY:4:1101:16532:1000 87 -997713494 3 150M -997713677 -333 AAAGTAACTTAATCTTTTAATCCTTTTTAGGATCATTAAACCAAATATTGCTGTAAAAAAAAAAGACAGTTATTCTATATTTTATCCTTGTCACCCCAACCAAATATTCTAATATAATACTTTTAAATAAAACTAACTACTACTTTATAA FFF:FFF::FFFFFF:FFFFFFFF:FF:FFFFFFFFFFFF:FFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFF:FFF:FFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFF NH:i:2 HI:i:1 AS:i:298 nM:i:0 A00709:52:H2KFFDSXY:4:1101:16532:1000 355 SEQ13 384256921 3 150M = 384257104 333 TTATAAAGTAGTAGTTAGTTTTATTTAAAAGTATTATATTAGAATATTTGGTTGGGGTGACAAGGATAAAATATAGAATAACTGTCTTTTTTTTTTACAGCAATATTTGGTTTAATGATCCTAAAAAGGATTAAAAGATTAAGTTACTTT FFFFFFFFFFFFFFFFFFFFFFFFFFFFF,FFFFFFFFF:FFF:FFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,FFFFFF:FFFFFFFFFFFF:FF:FFFFFFFF:FFFFFF::FFF:FFF NH:i:2 HI:i:2 AS:i:298 nM:i:0 A00709:52:H2KFFDSXY:4:1101:16532:1000 403 SEQ13 384257104 3 150M = 384256921 -333 TATCGACAAGAAAGTTTGCGACCTCGATGTTGAATTAAAAGTTCTTTATAGAGTAGAGACTATAATAGAAGGTCTGTTCGACCTTTAAAATTTTACATGATTTGAGTTCAGACCGGTGTAAGCCAGGTTGGTTTCTATCTTTCAGTAATT FFFFFFFFFFFF:F:,F,::FF:FFFFF:F:FFFFFF:FFFFFFFFFFFFFFFFFFFFF,:FFFF:FFFFFFF,FFFFFFFFFF:FFFFFF,FFFFFFFF:FFFF:F[W::sam_parse1] mapped query cannot have zero coordinate; treated as unmapped
Hope there is something that could be done. Best regards, Luka
alexdobin
commented
4 years ago Hi Luka,
thanks! If it appears in SAM, it looks like a bug... At the genome generation, your first post says you used --genomeChrBinNbits 898 Is this a typo? Was it --genomeChrBinNbits 8? 898 could have created the problems.
Cheers Alex
ljudevitluka
commented
4 years ago Hi Alex, I used --genomeChrBinNbits 898, I calculated it as genomesize/numberofscaffolds and that is around 898. Did I do something wrong there, should I use --genomeChrBinNbits 8 instead? If it helps here is a link to the webpage of genome that I am using: http://marmorkrebs.dkfz.de/
In the meantime I did everything again with --genomeChrBinNbits 12. And it seems that everything worked fine. The mapping results (% of mapped reads) are the some, the only problem is that the run takes 8 hours instead of 20 minutes, as you can see from progress file, is there a way to optimize that? M/hr number length unique length MMrate multi multi+ MM short other Apr 25 13:11:26 0.7 87466 295 74.1% 284.4 0.5% 15.5% 0.2% 0.0% 10.0% 0.2% Apr 25 13:18:32 3.3 787853 295 74.2% 284.3 0.5% 15.5% 0.2% 0.0% 10.0% 0.1% Apr 25 13:25:33 4.2 1486064 295 74.1% 284.4 0.5% 15.6% 0.2% 0.0% 10.0% 0.1% ..... Apr 25 21:03:57 5.8 46757697 294 74.1% 283.9 0.5% 15.6% 0.2% 0.0% 10.0% 0.2% Apr 25 21:04:57 5.9 47192970 294 74.1% 283.9 0.5% 15.6% 0.2% 0.0% 10.0% 0.2% Apr 25 21:08:06 5.9 47504246 294 74.1% 283.9 0.5% 15.6% 0.2% 0.0% 10.0% 0.2% ALL DONE!
Best regards and thank you for your suggestion it saved me,
Luka
alexdobin
commented
4 years ago Hi Luka,
--genomeChrBinNbits 12 should be fine To increase speed, you can try to reduce --seedPerWindowNmax from default 50 to 30 or even 20.
Cheers Alex
Hello,
I am working on RNA-seq data, that consists of 150 bp paired-end reads and I was hoping to align them to the already published reference genome of the some species. Genome size is: ~3300000000 bp, with more then 3750000 scaffolds/contigs and fragments ranging from 100 to 770000 bp, N50= 31000. So it's a fragmented genome, but only one that is available. I successfully (more or less) mapped my reads against it but unfortunately I can not make an index file (.bai) form the .bam output that I need to visualize the data in IGV.
Code I used to run STAR: STAR runMode genomeGenerate --runThreadN 8 --genomeDir Star_Genome --genomeFastaFiles genomefile.fa --genomeChrBinNbits 898
STAR --genomeDir Star_Genome --readFilesIn R21_forward_paired.fq.gz R21_reverse_paired.fq.gz --runThreadN 8 --alignIntronMax INTRONMAX --outSAMstrandField intronMotif --readFilesCommand zcat --outFileNamePrefix R21_STAR --outReadsUnmapped Fastx --outSAMtype BAM SortedByCoordinate
I tried samtools index and index -c, also changing -m values, it didn't work with the following error: [E::hts_idx_push] Region 536930388..536930896 cannot be stored in a bai index. Try using a csi index with min_shift = 14, n_lvls >= 6 samtools index: failed to create index for "NOGFFAligned.sortedByCoord.out.bam": Numerical result out of range Then I tried sambamba, got the index, but when I load it to the IGV I see no mapped reads, and I should see many of them. Samtools suggests that this is a problem with aligner, the one thing that confused me are the negative coordinates on multiple reads: A00709:52:H2KFFDSXY:4:2472:30074:12806 163 SEQ13 -2147377808 255 150M = -2147377704 158 GCAAGGCGAGGTTGCAGCTGCTTTTAAAAGGGGGGAGAAACCAAGTGATTATTTCAGTAGTGACGTGGGCTGTCAGAAGGAAGATTGGAGTCGGTGCCTCCAGCTCCTCGTAGCTCAGACGTCTTCCGGAGCCCCACGAGGCTGGTGATC FFFFFFFFFFF:FFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFF NH:i:1 HI:i:1 AS:i:202 nM:i:0
As for the mapping results, they look really good which makes thing even more confusing :/ Started job on | Apr 20 18:57:50 Started mapping on | Apr 20 18:58:08 Finished on | Apr 20 19:13:43 Mapping speed, Million of reads per hour | 182.90
Number of reads unmapped: too many mismatches | 0 % of reads unmapped: too many mismatches | 0.00% Number of reads unmapped: too short | 4742645 % of reads unmapped: too short | 9.98% Number of reads unmapped: other | 72637 % of reads unmapped: other | 0.15% CHIMERIC READS: Number of chimeric reads | 0 % of chimeric reads | 0.00%
Since I am new to bio informatics I am not sure what to do, and how to make this work and get the output to IGV. Also, I tried everything form scratch with GFF file (building the genome index mapping the reads and still I could not get the .bai file). Could you please tell me how to proceed?
Many thanks and best regards,
Luka