graphmap spliced alpha: no spliced alignments reported

[aechchiki]

Hi Ivan,

we tested Graphmap to align RNA-seq reads to the reference genome. For benchmark purposes, we tried to align the reference transcripts in fasta format on the reference genome, both from Ensembl.

I tried to install and run Graphmap as described in the docs for spliced alignments. I compiled the rna-alpha branch, as suggested, and run this command: /home/aechchik/graphmap/bin/graphmap-not_release align -x rnaseq -t 16 -r reference-genome.fa -d reference-transcripts.fasta -o graphmap_ref.sam

The main problems is that the local alignments for the transcripts sub-features are not reconciled: if a transcript has two exons, then both are reported but their local alignments are not merged to a spliced alignment, reported separately as multiple alignments instead. For example - transcript FBtr0344900 has two annotated exons:

cat Drosophila_melanogaster.BDGP6.84.gtf | grep FBtr0344900 | grep exon
4    FlyBase    exon    42774    43157    .    +    .    gene_id "FBgn0266617"; transcript_id "FBtr0344900"; exon_number "1"; gene_name "CR45124"; gene_source "FlyBase"; gene_biotype "lincRNA"; transcript_name "CR45124-RA"; transcript_source "FlyBase"; transcript_biotype "lincRNA"; exon_id "FBtr0344900-E1";
4    FlyBase    exon    43231    43374    .    +    .    gene_id "FBgn0266617"; transcript_id "FBtr0344900"; exon_number "2"; gene_name "CR45124"; gene_source "FlyBase"; gene_biotype "lincRNA"; transcript_name "CR45124-RA"; transcript_source "FlyBase"; transcript_biotype "lincRNA"; exon_id "FBtr0344900-E2";

In the alignment file by graphmap, it figures twice as local alignment - once as mapping on the forward strand, once as secondary mapping on the forward strand ($2):

cat graphmap_ref.sam | grep FBtr0344900
FBtr0344900 0   4   42774   40  384M144S    *   0   0TGCGACATTGTTCTACGATGACTACAAAAAATGACCAATAACTTCTATAAACCAATACGATATGTCAGGAGTTTCGGTCCCATACGAAGTCGCCGACTTAAGTATTTTATTTTTATTTTGATATGTGTTTGCTATTTTACCTTGTCGAATGCTTCCACACGCTATGAGAATACCATCGTGAGCGTAGCTTACTACTAGAATTTTGTTGAAGTTATTGACAAGCGATGTCTCAATATCTTCCGGACAGCCTCCAGCGTGACATTGCGGGGAATCATGTAACGGCCCAGTAACAGCCTCGGCCAGCACTCGAAGGTTTTCGTTAAGTTTAAGTATTTTATTTGTAGCACCCGCAAACAAAACATTGTGCATAAAGTCGAAGCTCATCTGGAAGCTGTTGATTGAACTGGTATTGATGGCAAGTTAAACTGGGCGACTATGTCATTTAAGGGAGATAACGCCTGAGCCGGCAGTTCTTCAATGCAGTTAACGCAATAATGCTGAGAACCGAGTATGATAATAATACACAGT   *   MD:Z:384    NM:i:0  AS:i:1920   H0:i:1  ZE:f:4.54542e-127   ZF:f:0.299417   ZQ:i:528ZR:i:1348131
FBtr0344900 256 4   43231   40  384S143M1S  *   0   0TGCGACATTGTTCTACGATGACTACAAAAAATGACCAATAACTTCTATAAACCAATACGATATGTCAGGAGTTTCGGTCCCATACGAAGTCGCCGACTTAAGTATTTTATTTTTATTTTGATATGTGTTTGCTATTTTACCTTGTCGAATGCTTCCACACGCTATGAGAATACCATCGTGAGCGTAGCTTACTACTAGAATTTTGTTGAAGTTATTGACAAGCGATGTCTCAATATCTTCCGGACAGCCTCCAGCGTGACATTGCGGGGAATCATGTAACGGCCCAGTAACAGCCTCGGCCAGCACTCGAAGGTTTTCGTTAAGTTTAAGTATTTTATTTGTAGCACCCGCAAACAAAACATTGTGCATAAAGTCGAAGCTCATCTGGAAGCTGTTGATTGAACTGGTATTGATGGCAAGTTAAACTGGGCGACTATGTCATTTAAGGGAGATAACGCCTGAGCCGGCAGTTCTTCAATGCAGTTAACGCAATAATGCTGAGAACCGAGTATGATAATAATACACAGT   *   MD:Z:143    NM:i:0  AS:i:715    H0:i:1  ZE:f:3.39258e-42    ZF:f:0.299417   ZQ:i:528ZR:i:1348131

What we would rather see instead would be something like the output of STAR:

cat star_ref.sam | grep FBtr0344900
FBtr0344900    0    4    42774    255    384M73N144M    *    0    0TGCGACATTGTTCTACGATGACTACAAAAAATGACCAATAACTTCTATAAACCAATACGATATGTCAGGAGTTTCGGTCCCATACGAAGTCGCCGACTTAAGTATTTTATTTTTATTTTGATATGTGTTTGCTATTTTACCTTGTCGAATGCTTCCACACGCTATGAGAATACCATCGTGAGCGTAGCTTACTACTAGAATTTTGTTGAAGTTATTGACAAGCGATGTCTCAATATCTTCCGGACAGCCTCCAGCGTGACATTGCGGGGAATCATGTAACGGCCCAGTAACAGCCTCGGCCAGCACTCGAAGGTTTTCGTTAAGTTTAAGTATTTTATTTGTAGCACCCGCAAACAAAACATTGTGCATAAAGTCGAAGCTCATCTGGAAGCTGTTGATTGAACTGGTATTGATGGCAAGTTAAACTGGGCGACTATGTCATTTAAGGGAGATAACGCCTGAGCCGGCAGTTCTTCAATGCAGTTAACGCAATAATGCTGAGAACCGAGTATGATAATAATACACAGT    *    NH:i:1    HI:i:1    NM:i:0    MD:Z:528

Am I doing something wrong or simply this feature is not implemented yet?

Another strange behaviour as seen in the output alignment file is that some of the reads ID appear once in the alignment file, but there is no other alignment of the same read ID reporting the primary alignment flag. For example: cat graphmap_ref.sam | grep -v '^@' | awk '$2==256 {print $0}' | head -1 | cut -f1 outputs: FBtr0300689. If we grep this ID out of the mapping file, then we are returned one hit only:

cat graphmap_ref.sam | grep -v '^@' | grep FBtr0300689
FBtr0300689 256 2L  8191    40  586S1293M1S *   0   0CTACTCGCATGTAGAGATTTCCACTTATGTTTTCTCTACTTTCAGCAACCGAGAAGAGAACCCACGTTTGAACAAGTATCGGCGTGTGGACAACAGCTATCCCCGCTTCATAACGAATGAGGCTGCCGAGGACCTGATTTACAAGAAGTCCATGGGCGAGCGGGATCAGCCACAGAGCTCAGAGCGGATCTCAATATTTAATCCGCCAGTATACACGCAGCACCAGGTGCGCAATGAAGCCCCCTACATACCCACCACATTTGACCTCCTCTCAGACGATGAGGAGTCGTCACAGAGAGTTGCCAACGCCGGGCCATCTTTCAGGCCCTTGACTTACTCGGATGCTGTGCGTCTAAGCCAGAATGGCTTCGCCAACTCCCGCGTAAGTGGGCACTCCAGCTATACGGTGCGCAGACCACCGGCACTAGTTGACAGAAGCATTCTATCCCAGGAAATGGAGCGCATGGACCAAGAGCAGTATATCTACCTTATCCGTACCGCAGCCCAAAGTAATTCCGTGGGCAGTCACTACGCCGAACCGGTTACTGATAACTCGGAGGTCAAGAAAGTCAGTGAAACCAACAAAAGTGACGCACCACAGCCGTTAACCCCTCAACCTACCAGACTCACCAGAACAGAATCCTTGCACCGTCGTTTTGCCAGCTGCGTCAACTTAAATGATGACTTCGCCGAGCAATTTAAAGCAAGAGCGGCGGACTGTGAAGAGAAATCCAAACATCGTCTTAGATTAGCTGAAGAGCAGAGGCTTTTTTCGAATTTCAGTGCTATAAAGAACATAGATGAACTCCGTGCCTATGAACGAAAAGTAGTGGAAAACATATTCCAGTCTTGTATCGCCCACAAGCCCATTTTTGTACTCGGGCCCTTGGACAAGCCAAATGTGAAGAAAGTGACCAAGCTCATTCCGTTAACAGAGGAGCACCACGATCGCTTTAACGAAATTACACAGGATGATAAATCGACGGTATGGCAACGAATATATTGATGTCTTTCGTACCCATTGAAAACGTTGTGGTGCTTGCGCTTTAAAATCTTATATTAGGAAATTATTTTTAAATTTAACCTACACATAACTACCGAAGACATATGCACGTTTATTAATGGGAAATGGCTTAACGACGAGGTCATTAACTTTTACATGTCCTTGCTGACAGAACGGTCGGAGAAGAGATCTGGCGTACTTCCCGCCACTTACGCCATAAACACATTCTTCGTGCCCCGCCTCCTGCAAGCTGGGCATGCAGGCATTAAGCGCTGGACTCGCAAAGTGGACTTGTTCAGCAAGGACATAATCCCGGTACCAGTGCACTGCAACGGCGTCCACTGGTGCATGGCCATCATACACTTGCGGAACAAGACAATCCGGTATTATGACTCAAAGGGAAAGCCAAACCGACCAGTGCTGGACGCTCTAGAGAAATATCTACGCGAAGAGTCAATATTCAAGCCCAAAAAGCAGTTTGATACCAGCGATTTTGTTATTGAGAGCGTGCAGAATATACCACGACAGTTAGATGGCAGCGATTGCGGTATCTTCAGCTGCATGTTCGCCGAGTATATAACGTGTGATGTGCCAATTACCTTTACCCAGTCGGAAATGTTGTACTTCCGCAAGAAGATGGCTCTAGAAATCGTCGACGGAGAGTTGTGACAGTAGAATCACACAGCTACGCAAGAATGTGGAGAATCCAGTTTAGTTATTTTTACAAATCTTACGTAAACACTCCAAGCATGAATTCGCAACAAGTGCTTAGCTATTTAATTGAATTGAGCTGGCCGAGAGATGTGCTGGTGCAATAACTTGTTCTCATATCTGATTGTAACAGAGAATCTAGTTTTTCAATAAAATTTCCCCAAGT   *   MD:Z:1293   NM:i:0  AS:i:6465   H0:i:1  ZE:f:0  ZF:f:0.276046   ZQ:i:1880   ZR:i:23513712

In my opinion, this read ID should have alignment flag equal to 0 instead.

I compiled the branch alpha on CentOS 6.9:

aechchik@frt:~/software/graphmap$ git branch
  master
* rna-alpha

And graphmap-not_release is the only executable in the execution folder:

aechchik@frt:~/software/graphmap/bin$ ls
graphmap-not_release

Cheers, Amina

[isovic]

Dear Amina,

Thank you for testing and reporting back! You are right about exons being reported as separate SAM alignments. This was an intentional approach for the proof of concept. However, we are working on implementing a solution which would produce alignments similar to STAR, but it might take some time to finalize this.

Regarding the missing primary alignments - that should definitely not happen. Could you by any chance share a sample of your data which reproduces this issue, so I can take a look?

Best regards, Ivan.

[aechchiki]

Hi Ivan,

I ran the following command: /home/aechchik/graphmap/bin/graphmap-not_release align -x rnaseq -t 16 -r reference-genome.fa -d reference-transcripts.fasta -o graphmap_ref.sam where:

• reference-genome.fa is a concatenation of chromosomal genome files at: ftp://ftp.ensembl.org/pub/release-89/fasta/drosophila_melanogaster/dna/ : cat Drosophila_melanogaster.BDGP6.dna.chromosome.*.fa.gz
• reference-transcripts.fasta is an extraction of the fasta sequences from the reference annotation at: ftp://ftp.ensembl.org/pub/release-89/gtf/drosophila_melanogaster/Drosophila_melanogaster.BDGP6.89.gtf.gz using the gffread utility, using the command: gffread -w reference_transcripts.fasta -g reference_genome.fa reference_transcripts.gtf You should be able to reproduce exactly what I got.

(I edited the comment, for some reason the links to the Ensembl ftp were not displayed )