cufflinks, transcripts extraction

[kenietz]

Hi guys,
some clients want me to extract the consensus sequence of the transcripts assembled with cufflinks. Well i've searched around and found a way how to extract the consensus using mpileup. But my problem now is how to know which transcript has some consensus at all.

For example i run cufflinks on a bam file and i get some tracking files and transcripts.gtf. So far so good. But then i open the genes_fpkm.tracking and choose for example the following genes:

gene4 - - gene4 Gm10568 - chr1:3680154-3681788 - - 0 0 0 OK

gene9 - - gene9 Gm7369 - chr1:4591096-4611406 - - 0 0 0 OK

I take the coordinates and for the first one i get no consensus sequence but for the second i get something. When i go to IGV with those coordinates i see no reads for the first gene and one read for the second.

So what shoud i do? Just skip those lines with zeros or not?

Btw, i just remembered another thing. With this one liner provided by samtools man page for extracting a consensus i get all the Ns in the beginning of the chromosome. I mean with that line:

samtools mpileup -r REGION -uf REFGENOME bam | bcftools view -cg -| vcfutils.pl vcf2fq > cns.fq

Someone knows how to get rid of it without using additional programming?

For now i am using for the extraction a small perl script parsing the vcf file coming out of these 3 lines:

$res=`samtools mpileup -r $region -uf $refgen $bam_in > $bcf1`;
$res=`bcftools view -bcg $bcf1 $region > $bcf2`;
$res=`bcftools view $bcf2 > $vcf`;

Thanks for any tips and if someone is interested in the script itself please say
Cheers