This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Visualizing alignment of kallisto .bam files to hg19 (and custom plasmids)

Hello,

I'm trying to visualize kallisto read alignment to a human genome plus an additional sequence using IGV. I can produce .bam and .bai files using kallisto, but when I open them in IGV alongside the built-in "hg19" option, I don't see anything:

example IGV window - no alignment where I would expect to see something

Gathering Files

First, I downloaded my reference human genome from here, along with a gtf:

http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/

wget 'http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/mrna.fa.gz'

Gtf downloaded from the UCSC genome browser using these instructions: hg38 annotation file (gtf) UCSC format

Then, I download my additional sequence, "plas_shp", which looks like this:

head plasmidFasta/plas_shp_up.fa
>plas_shp
GTCGACATTGATTATTGACTAGTTATTAATAGTAATCAATTACGGGGTCATTAGTTCATA
GCCCATATATGGAGTTCCGCGTTACATAACTTACGGTAAATGGCCCGCCTGGCTGACCGC
CCAACGACCCCCGCCCATTGACGTCAATAATGACGTATGTTCCCATAGTAACGCCAATAG
GGACTTTCCATTGACGTCAATGGGTGGAGTATTTACGGTAAACTGCCCACTTGGCAGTAC
ATCAAGTGTATCATATGCCAAGTACGCCCCCTATTGACGTCAATGACGGTAAATGGCCCG
CCTGGCATTATGCCCAGTACATGACCTTATGGGACTTTCCTACTTGGCAGTACATCTACG
TATTAGTCATCGCTATTACCATGGTCGAGGTGAGCCCCACGTTCTGCTTCACTCTCCCCA
TCTCCCCCCCCTCCCCACCCCCAATTTTGTATTTATTTATTTTTTAATTATTTTGTGCAG
CGATGGGGGCGGGGGGGGGGGGGGGGCGCGCGCCAGGCGGGGCGGGGCGGGGCGAGGGGC

Producing Index

I begin by producing my kallisto index. The code looks like this:

#declare Variables                                                                                                                          
humFa="myPath/hg19_transcriptome/*.fa"  
plasFa="myPath/plasmidFasta/*.fa"

 # run index gen
  kallisto index -i /myPath/humPlas_kallisto_transcripts_plasmids_hg19.idx $humFa $plasFa --make-unique

Producing .bam file

I modified the gtf file I downloaded above to add my extra sequences, so it looked like this:

tail myGtf.gtf
chr19_gl000209_random   hg19_refGene    stop_codon  145118  145120  0.000000    +   .   gene_id "NM_012312"; transcript_id "NM_012312"; 
chr19_gl000209_random   hg19_refGene    exon    145079  145743  0.000000    +   .   gene_id "NM_012312"; transcript_id "NM_012312"; 
plas_shp    AddedGenes  exon    1   12915   .   +   0   gene_id "plas_shp"; gene_name "plas_shp"; transcript_id "plas_shp"; tss_id "plas_shp";

Then I aligned kallisto to the index I produced above using the above:

# run index gen
kallisto quant -i $kallistoIdx -o $outputFileLoc/$sample --genomebam --gtf $myGtf $Read1 $Read2 | samtools view -bS > $sample_align.bam

When I convert the output .bam file to .sam, it looks like this:

NS500127:27:HW5JKBGXY:1:11101:17213:1141    131 *   0   0   76M *   0   0   GGCGTTGCTGCGGGTGGCCTGCGTGTGCCGCTTATGGAGGGAGTGTGTGCGCAGAGTATTGCGGACCCATCGGAGC    AAAAAEEEEEAEEEAEEEEE/AEE/AEEEEEEAEEEE/EEEEE/EAE/AAE//E/AEE/EEAEA//E/<A/EEEEA    ZW:f:0
NS500127:27:HW5JKBGXY:1:11101:24672:1142    67  *   0   0   76M *   0   0   TTTGATCTTCTTAGTTTTCTTCTTCTTATCCTTACCGCTGTCATCCTCCTCATCTGAACCCACATCTTCGATCTTG    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAEEEEEEEAEEEEEEEEEAE/EEEEEEEEEEEE//EEEEA    ZW:f:0

The problem

When I try to open the .bam file in IGV, I don't see any alignments. I've tried doing the above with hg19/GRCh37 and 38 files from different locations. I've tried including and excluding the chromSizes file. I'm not sure what else to try.

Is there anything immediately obviously problematic about my approach?

Is there a viewer that might be more forgiving than IGV when it comes to visualizing kallisto .bam files?

Has anyone successfully managed to visualize their kallisto .bam files in IGV before?

Thank you for your help!

igv kallisto rna-seq
• 5,252 views • link
updated by pmelsted pmelsted • written by Kristin Muench Kristin Muench

2 answers

Everything up to the GTF file was correct. The only thing missing is that kallisto requires a gene entry for each transcript and the exons of each transcript. Also it's important the the field separators for the GTF are TAB characters, where as the key-value entries are separted by space and ;

plas_shp    AddedGenes  gene    1   12915   .   +   0   gene_id "plas_shp"; gene_name "plas_shp";
plas_shp    AddedGenes  transcript  1   12915   .   +   0   gene_id "plas_shp"; gene_name "plas_shp"; transcript_id "plas_shp";
plas_shp    AddedGenes  exon    1   12915   .   +   0   gene_id "plas_shp"; gene_name "plas_shp"; transcript_id "plas_shp"; tss_id "plas_shp";

this line kallisto quant -i $kallistoIdx -o $outputFileLoc/$sample --genomebam --gtf $myGtf $Read1 $Read2 | samtools view -bS > $sample_align.bam is unnecessary since kallisto with --genomebam doesn't write to standard output, but produces a sorted indexed BAM file under outputFolder/pseudoalignments.bam.

After this you should be able to view this using IGV. You need to either index your plasmid fasta file using samtools faidx or you can add them as extra chromosomes to the Human reference (but give them a new chromosome name and modify the GTF and run samtools faidx again).

I tested this on a small example 

@r1
TGACTAGTTATTAATAGTAATCAATTACGGGGTCATTAGTTCATAGCCCATATATGGAGT
+
AAAAAEEEEEAEEEAEEEEE/AEE/AEEEEEEAEEEE/EEEEE/EAE/AAE//E/AEE/E
@r2
TGGGCATAATGCCAGGCGGGCCATTTACCGTCATTGACGTCAATA
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAEEEEEEEA

using the head of the plasmid sequence you provided and got the following results kallisto quant -i ind -o out --genomebam --gtf plas.gtf --single -l 200 -s20 r1.fastq and looking at the output samtools view out/pseudoalignments.bam

r1  0   plas_shp    16  255 60M *   0   0   TGACTAGTTATTAATAGTAATCAATTACGGGGTCATTAGTTCATAGCCCATATATGGAGT    AAAAAEEEEEAEEEAEEEEE/AEE/AEEEEEEAEEEE/EEEEE/EAE/AAE//E/AEE/E    ZW:f:1
r2  16  plas_shp    272 255 45M *   0   0   TATTGACGTCAATGACGGTAAATGGCCCGCCTGGCATTATGCCCA   AEEEEEEEAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA   ZW:f:1
• 1 views • link
written by pmelsted pmelsted

Thank you for this extremely helpful answer! I'll give it a try!

• 1 views • link
written by Kristin Muench Kristin Muench

Thank you again for your suggested tips. I tried editing my gtf, removing the samview portion of my kallisto quant line, and re-running the pipeline. This time, I did it both the original way (trying to build index from plasmid+ transcriptome) as well as a simplified vdrsion (trying to build index from transcriptome .fa only). Unfortunately, both versions failed to align properly (e.g. sam output very similar to what I posted above with no chromosome names) and I could not visualize them on IGV. I'll focus on figuring out what's wrong with just the basic transcriptome alignment.

Regarding your other suggestion using samtools faidx - could you help me understand how this would enter into the pipeline? I see that it makes .fa.fai files, but if you aren't feeding .fai files into kallisto index, just .fa files (see section "producing index" for example), why would it matter they exist?

Thank you again for your help!

• 1 views • link
written by Kristin Muench Kristin Muench

Your bam has no indication whatsoever what chromosome each read goes to. So of course IGV can't put reads on the human chromosomes.

If you want to look at alignments in IGV, you need to align to the human genome with a real aligner, like STAR, not just count transcripts with kallisto.

• 1 views • link
written by swbarnes2 swbarnes2

Thanks for pointing out where it's falling down. You actually can use IGV with kallisto now: https://pachterlab.github.io/kallisto/pseudobam.html

• 1 views • link
written by Kristin Muench Kristin Muench

Thanks for pointing out where it's falling down. You actually can use IGV with kallisto now: https://pachterlab.github.io/kallisto/pseudobam.html

• 1 views • link
written by Kristin Muench Kristin Muench

Thanks for pointing out where it's falling down. You actually can use IGV with kallisto now: https://pachterlab.github.io/kallisto/pseudobam.html

• 1 views • link
written by Kristin Muench Kristin Muench

Log in to answer this question.

More posts like this
  • feature count command error
    written by Adyasha • Adyasha
    Hi everyone , I have ran this feature count command : featureCounts -T 4 -p -a gencode.v43.basic.annotation.gtf -o featurecount.txt *.bam this gave me this error …
  • kallisto genomebam error
    written by Stella • Stella
    Hello! I am trying to produce bam files to load to igv after kallisto quant with --genobam option. This is the code that I am …
  • Should kallisto take hours to run?
    written by Kristin Muench Kristin Muench
    Hello, I'm trying to align some test .fastq files to a reference using kallisto. The reference is composed of human transcriptome + a couple of …
  • Adding plasmid sequences to gtf file
    written by Kristin Muench Kristin Muench
    Hello, I am trying to align some .fastq files to both the hg19 human genome build AND a couple of plasmid sequences. I am trying …
  • Ht-seq count for ncRNA.fasta. count=0
    written by call • call
    Hello. I used a ncRNA.fa to alignment using RNA-seq samples in TOPHAT *FASTA: >ENSSSCT00000024762.1 ncrna:novel chromosome:Sscrofa10.2:2:114672074:114672179:1 gene:ENSSSCG00000023111.1 gene_biotype:snRNA transcript_biotype:snRNA gene_symbol:U6 description:U6 spliceosomal RNA [Source:RFAM;Acc:RF00026] GCGTCACTGTAGTAGCACACAGACTAAAACTGGAATGATACAGAGAAGATTAGCATAGAT …
  • Htseq-Count Error With Strand-Specific Data
    written by Nicolas Rosewick Nicolas Rosewick
    <p>Hi,</p> <p>I have 2x50 bp human strand-specific data. I aligned them against hg19 with STAR and after that I extracted the read count per gene …