Visualizing Alignments to a gene

[mgibson]

I will apologize in advance if this is a really obvious question. I am analyzing some RNA seq data using TopHat/Bowtie for alignment and Cufflinks for expression analysis. I want to pull out the consensus sequence for all the reads that align to a particular gene. I know that I can use SAMtools to pull out all the raw reads that align to a region, but that is still thousands of individual reads to deal with.

Is there a function that I am missing in the software packages I currently use? Or is there some alternative piece of software that I should be using?

Thanks in advance!

[gringer]

Galaxy has a nice pileup-to-interval function that takes as input the mpileup output of samtools, and generates a consensus sequence for regions above a certain coverage. Tablet is nice for visualising mapping and the consensus sequence and exon boundaries:

http://bioinf.scri.ac.uk/tablet/

[Richard Finney]

If you just want to view chunks of RNA seq bams, try this script, it creates a wig (wiggle) file which you can upload to UCSC browser as a custom track. In this case, it is called "bam2wig"

#example usage: ./bam2wig example.bam chr5:10333595-15399999 fileName.wig
#this uses SAMT variable, you must edit export line (bash specific ) to point to your samtools executable

export SAMT=samtools-0.1.18/samtools
echo "track type=wiggle_0 name=fileName.wig description=examplecustomwig" > $3
$SAMT depth -r $2 $1 | awk '{if ($1!=p) {print "variableStep chrom="$1;} print $2" "$3;p=$1;}' >> $3

[dglemay]

Thanks

Dear Richard,
Great script. Thanks for sharing!

For others who use it, don't forget to index your bam file first
samtools index <yourfile>.bam

-Danielle

[jflowers]

If you want to stick with samtools, you can use bcftools and vcfutils.pl vcf2fq (which come with your samtools installation) to get a consensus. This can all be done in one command line, like:

samtools view -u aln.bam chromsome_name:X1-X2 | samtools mpileup -uf ref.fa - | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq

Where aln.bam is your alignment, chromosome_name is the identifier in your original ref.fa, and X1 and X2 are the coordinates of interest.

See documentation here:
http://samtools.sourceforge.net/mpileup.shtml