SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
strand-specific cDNA papori General 26 08-02-2017 11:30 AM
Strand-specific library appears not strand-specific oligo Illumina/Solexa 7 12-08-2011 10:54 AM
strand specific RNAseq Pepe Sample Prep / Library Generation 5 10-27-2011 08:58 AM
cuffdiff not strand-specific? burkard Bioinformatics 0 11-15-2010 11:30 AM
strand-specific Illumina protocol PFS Bioinformatics 1 06-28-2010 08:53 AM

Reply
 
Thread Tools
Old 07-27-2011, 11:47 PM   #1
frymor
Senior Member
 
Location: Germany

Join Date: May 2010
Posts: 149
Unhappy strand specific wig files

Hi,

I am trying to create strand-specific wig files from a bam file.
I used the pileup option and i thought it worked very good. Now I have checked it again with the mpileup and it looks very strange to me.

This is how I did it before:
Code:
samtools view -h -f 0x10 dilpAll.sorted.bam | samtools pileup -S - | cut -f 1,2,4| awk 'BEGIN{print "tracktype=\"wiggle_0\"\tname=\"dilp (reversestrand)\"\tdescription=\"dilp (reverse strand)\""} {if(NR==1){curChr=$1; print "variableStep chrom="curChr" span=1";} 
if($1 == curChr){print $2"\t"$3;}else{curChr=$1; print "curChr="curChr;}}' > wigFiles/dilp_rs.wig

samtools view -h -F 0x10 dilpAll.sorted.bam | samtools pileup -S - | cut -f 1,2,4| awk 'BEGIN{print "track type=\"wiggle_0\"\tname=\"dilp (forward strand)\"\tdescription=\"dilp (forward strand)\""} {if(NR==1){curChr=$1; print "variableStep chrom="curChr" span=1";
} if($1 == curChr){print $2"\t"$3;}else{curChr=$1; print "curChr="curChr;}}' > wigFiles/dilp_fs.wig
Now I just have changed it a bit to fit the mpileup options:
Code:
samtools view -b -f 0x10 dilpAll.sorted.bam | samtools mpileup  - | cut -f 1,2,4| awk 'BEGIN{print "tracktype=\"wiggle_0\"\tname=\"dilp (reversestrand)\"\tdescription=\"dilp (reverse strand)\""} {if(NR==1){curChr=$1; print "variableStep chrom="curChr" span=1";} 
if($1 == curChr){print $2"\t"$3;}else{curChr=$1; print "curChr="curChr;}}' > wigFiles/dilp_rs.wig

samtools view -b -F 0x10 dilpAll.sorted.bam | samtools mpileup - | cut -f 1,2,4| awk 'BEGIN{print "track type=\"wiggle_0\"\tname=\"dilp (forward strand)\"\tdescription=\"dilp (forward strand)\""} {if(NR==1){curChr=$1; print "variableStep chrom="curChr" span=1";
} if($1 == curChr){print $2"\t"$3;}else{curChr=$1; print "curChr="curChr;}}' > wigFiles/dilp_fs.wig
Than I compared it with a wig file for the omplete bam file:
Code:
samtools mpileup dilpMerged.sorted.bam | perl -ne 'BEGIN{print "track type=wiggle_0 name=fileName description=fileName\n"};($c, $start, undef, $depth) = split; if ($c ne $lastC) {print "variableStep chrom=chr$c\n"; };$lastC=$c;next unless $. % 10 ==0;print "$star
t\t$depth\n" unless $depth<3;'  > wigFiles/dilpMerged.wig
When doing it for both strands I get more line than I have together in the wig file for the complete bam file.
-f: 113905
-F: 121720
w.o.: 141199

Does someone have an explanation for these numbers?

Thanks
Assa
frymor is offline   Reply With Quote
Reply

Tags
bam, convert, strand-specific, wig

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 01:59 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO