cufflinks, transcripts extraction
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
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