Changes between Version 1 and Version 2 of SOPs/coordinates

Timestamp:

12/18/14 15:46:57 (11 years ago)

Author:

byuan

Comment:

--

SOPs/coordinates

@@ -151 +151 @@
 }
  }}}
+
 === Convert blat to gff ===
 
   * Use /nfs/BaRC_Public/BaRC_code/Perl/blat2gff/blat2gff.pl
+
   
  {{{
    Convert BLAT output file (PSL format) into GFF format (v1.1 14 Dec 2010)
     blat2gff.pl blatFile dataSource(ex:WIBR) > gffFile
-}}
+}}}
+
+
+=== create wiggle files for visualizing paired-end data mapping to the + and - strands ===
+
+1. split by strand by matched strand
+
+{{{
+# input:        accepted_hits.bam
+# output:       accepted_hits_negStrand.bam: mapped to negative strand
+#               accepted_hits_posStrand.bam: mapped to positive strand
+
+bsub "samtools view -f 16 -b accepted_hits.bam >| accepted_hits_negStrand.bam"
+bsub "samtools view -F 16 -b accepted_hits.bam >| accepted_hits_posStrand.bam"
+
+}}}
+
+2. split reads by pair
+{{{
+# input:        accepted_hits_posStrand.bam or accepted_hits_negStrand.bam
+# output:       1st pair: *_1stPair.bam
+#                       2nd pair: *_2ndPair.bam
+bsub "samtools view -b -f 0x0040 accepted_hits_posStrand.bam > accepted_hits_posStrand_1stPair.bam"
+bsub "samtools view -b -F 0x0040 accepted_hits_posStrand.bam > accepted_hits_posStrand_2ndPair.bam"
+bsub "samtools view -b -f 0x0040 accepted_hits_negStrand.bam > accepted_hits_negStrand_1stPair.bam"
+bsub "samtools view -b -F 0x0040 accepted_hits_negStrand.bam > accepted_hits_negStrand_2ndPair.bam"
+}}}
+
+3. convert from bam to bedgraph format
+{{{
+# input:        bam format: accepted_hits_*Strand_*Pair.bam
+#                       /nfs/genomes/mouse_gp_jul_07/anno/mm9.size: length of each chromosome, format like 
+#                                   chr1    197195432
+# output:       bedgraph format: accepted_hits_*Strand_*Pair.bedgraph
+bsub "genomeCoverageBed -split -bg -ibam accepted_hits_posStrand_1stPair.bam -g mm9.size >| accepted_hits_posStrand_1stPair.bedgraph"
+bsub "genomeCoverageBed -split -bg -ibam accepted_hits_posStrand_2ndPair.bam -g mm9.size >| accepted_hits_posStrand_2ndPair.bedgraph"
+bsub "genomeCoverageBed -split -bg -ibam accepted_hits_negStrand_1stPair.bam -g mm9.size >| accepted_hits_negStrand_1stPair.bedgraph"
+bsub "genomeCoverageBed -split -bg -ibam accepted_hits_negStrand_2ndPair.bam -g mm9.size >| accepted_hits_negStrand_2ndPair.bedgraph"
+}}}
+
+4. join the reads sharing the same strand
+{{{
+# This step is for fr-firststrand library (such as dUTP). which is
+    1+-,1-+,2++,2--
+
+    read1 mapped to ‘+’ strand indicates parental gene on ‘-‘ strand
+    read1 mapped to ‘-‘ strand indicates parental gene on ‘+’ strand
+    read2 mapped to ‘+’ strand indicates parental gene on ‘+’ strand
+    read2 mapped to ‘-‘ strand indicates parental gene on ‘-‘ strand
+ 
+# input:        bedgraph file from the same strand
+# output:       merged bedgraph: pos.bedgraph or neg.bedgraph
+unionBedGraphs -i accepted_hits_posStrand_2ndPair.bedgraph accepted_hits_negStrand_1stPair.bedgraph |awk '{ print $1"\t"$2"\t"$3"\t"$4+$5 }' >|pos.bedgraph
+unionBedGraphs -i accepted_hits_posStrand_1stPair.bedgraph accepted_hits_negStrand_2ndPair.bedgraph |awk '{ print $1"\t"$2"\t"$3"\t-"$4+$5 }' >|neg.bedgraph
+}}}
+
+  * To check how reads were stranded for strand-specific RNA-seq data, you can use infer_experiment.py from  [http://rseqc.sourceforge.net/ RSeQC]
+
+{{{
+infer_experiment.py -r mm9.refseq.bed12 -i accepted_hits.bam
+}}}
+
+5. convert bedgraph to bigwig
+
+{{{
+# get rid of header lines of mm9.size: the header line with "chrom   size" is removed
+# input:        mm9.size: length of each chromosome
+# output:       mm9.size_noHeader
+tail --line=+2 mm9.size > mm9.size_noHeader
+# convert bedgraph to bigwig
+# input:        bedgraph file: neg.bedgraph or pos.bedgraph
+#                       mm9.size_noHeader: length of each chromosome
+# *output:      bigwig format: neg.bw or pos.bw
+#                       neg.bw or pos.bw can be visualized with IGV/UCSC genome browser
+bsub bedGraphToBigWig neg.bedgraph mm9.size_noHeader neg.bw
+bsub bedGraphToBigWig pos.bedgraph mm9.size_noHeader pos.bw
+}}}