Question: Visualization for DE transcripts
0
4 months ago by
rbenel0
rbenel0 wrote:

I am trying to find a way to visualize the different transcripts that I found to be DE using this workflow (https://bioconductor.org/packages/release/workflows/vignettes/rnaseqDTU/inst/doc/rnaseqDTU.html).

I tried to re-run salmon and use the --writeMapping flag to receive a .sam file and then work from there, what is the best way to visualize a .sam file for pseudoaligned transcripts?

Thank You!

rnaseqdtu • 181 views
modified 4 months ago • written 4 months ago by rbenel0
0
4 months ago by
Michael Love26k
United States
Michael Love26k wrote:

I'll note that many of my UNC collaborators use STAR with --quantMode TranscriptomeSAM followed by Salmon for quantification. And this way you can view the genomic alignments, including those that don't map to the annotated transcriptome.

With Rob's group and Charlotte and Carl, we've worked recently on showing how for some reads, the genomics alignments are also helpful, although a combination appears the best option.

https://www.biorxiv.org/content/10.1101/657874v1

0
4 months ago by
rbenel0
rbenel0 wrote:

So in theory is isn't an issue to visualize the transcripts, if one follows a few steps.

1) Following using the --writeMapping flag when running salmon the pseudosam file needs to be sorted using samtools samtools sort -o /path/to/output.bam /path/to/input.sam

2) After sorting, the sorted .bam file needs to be indexed. samtools index -b path/to/sorted.bam

3) Finally, the transciptome needs to be indexed as well. samtools faidx path/to/gencode.v28.transcripts.fa -o path/to/gencode.v28.transcripts.fai

If when running salmon the --gencode flag was used, you will find that the gencode file that was used for indexing, when uploaded to IGV will not match the alignment file and this is because "the flag splits the transcript name at the first '|' character".

A quick and "dirty" fix for that is to do the same to the gencode fasta file, as such. awk 'BEGIN{FS="|"}{if(/^>/){print $1}else{print$0}}' ./gencode.v28.transcripts.fa > ./TranscriptNameGencode.v28.transcripts.fa

I hope this helps, and thanks to @rob-p @k3yavi @Michael Love and @daviesrob for answering questions along the way :)

1