Dinh 2011/NOTES/2011-11-14

From ZhangLabWiki
Jump to navigation Jump to search

Yi Zhang UNC-5caC/5fC[edit]

8,506,009 reads for 5fC 
39,354,845 reads for 5caC
12,518,460 reads for IgG
25,206,212 reads for Input

Mapped with bowtie: (# alignments)[edit]

bowtie -n 2 --best --strata -l 28 -k 1 -m 1 -p 16 --phred64-quals --un $f.unmappable.txt /media/1TB_store1/Rui_Scratch/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome $f $f.bowtie.map 2> $f.log 

> YZ_UNC-5caC-DIP.fastq.log <==
# reads with at least one reported alignment: 20985742 (53.32%)
# reads that failed to align: 2218974 (5.64%)
# reads with alignments suppressed due to -m: 16150129 (41.04%)
Reported 20985742 alignments to 1 output stream(s)
> YZ_UNC-5fC-DIP.fastq.log <==
# reads with at least one reported alignment: 3740642 (43.98%)
# reads that failed to align: 970771 (11.41%)
# reads with alignments suppressed due to -m: 3794596 (44.61%)
Reported 3740642 alignments to 1 output stream(s)
> YZ_UNC-IgG-DIP.fastq.log <==
# reads with at least one reported alignment: 8023681 (64.09%)
# reads that failed to align: 1072660 (8.57%)
# reads with alignments suppressed due to -m: 3422119 (27.34%)
Reported 8023681 alignments to 1 output stream(s)
> YZ_UNC-Input.fastq.log <==
# reads with at least one reported alignment: 16508323 (65.49%)
# reads that failed to align: 3733811 (14.81%)
# reads with alignments suppressed due to -m: 4964078 (19.69%)
Reported 16508323 alignments to 1 output stream(s)

Balance hits[edit]

2191  head -3740642 YZ_UNC-5caC-DIP.fastq.bowtie.map > sampled_YZ_UNC-5caC-DIP.fastq.bowtie.map &
2192  head -3740642 YZ_UNC-IgG-DIP.fastq.bowtie.map > sampled_YZ_UNC-IgG-DIP.fastq.bowtie.map &
2193  head -3740642 YZ_UNC-Input.fastq.bowtie.map > sampled_YZ_UNC-Input-DIP.fastq.bowtie.map &

Peak calling using MACS[edit]

macs14 -t sampled_YZ_UNC-5caC-DIP.fastq.bowtie.map -c sampled_YZ_UNC-IgG-DIP.fastq.bowtie.map -f BOWTIE -g mm -n 5caC-vsIgG -w --call-subpeaks > macs14_5caCvsIgG_balanced &
macs14 -t sampled_YZ_UNC-5caC-DIP.fastq.bowtie.map -c sampled_YZ_UNC-Input-DIP.fastq.bowtie.map -f BOWTIE -g mm -n 5caC-vsInput -w --call-subpeaks > macs14_5caCvsInput_balanced &

macs14 -t YZ_UNC-5fC-DIP.fastq.bowtie.map -c sampled_YZ_UNC-IgG-DIP.fastq.bowtie.map -f BOWTIE -g mm -n 5fC-vsIgG -w --call-subpeaks > macs14_5fCvsIgG_balanced &
macs14 -t YZ_UNC-5fC-DIP.fastq.bowtie.map -c sampled_YZ_UNC-Input-DIP.fastq.bowtie.map -f BOWTIE -g mm -n 5fC-vsInput -w --call-subpeaks > macs14_5fCvsInput_balanced &

Remapping using BWA to check enrichment in repeats[edit]

  • For ChIP-seq experiments, it is okay to use reads mapped to multiple locations, since a multi-mapped read can be split to its multiple locations, and we have data from the control experiment to normalize the enrichment results.
  • I decided to re-map using BWA since it have better sensitivity and I think better sensitivity is good for mapping to repetitive regions of the genome.
    • Use reads which map uniquely in one experiment -> unique genome test
    • Use reads which map up to many different regions in the genome -> repetitive genome test
  • Building mm9 bwa index: I copied genome.fa from Rui's iGenome folder:
cp /media/1TB_store1/Rui_Scratch/iGenome/Mus_musculus/UCSC/mm9/Sequence/WholeGenomeFasta/genome.fa BWAIndex/
/usr/local/bwa-0.5.9/bwa index -p mm9 -a bwtsw genome.fa &
  • Understanding BWA tags:
    • By default, aln finds all best hits and all suboptimal hits (within the given constraints)
    • Using samse, each read will be aligned and a best hit will be reported at random with the XA tag which will list up to a maximal number of alternative hits (default 3). Reads which map to more than (3) optimal+suboptimal locations will not have the XA tag but a best hit will be reported at random.
    • XT tag will indicate whether a read mapped uniquely or in repeats.
    • X0 tag will indicate how many best hits were found.
  • Understanding BWA flags:
awk '{print $2;}' *.sam | sort -u
    • 0 = mapped, hit given on forward strand
    • 4 = unmapped read
    • 16 = mapped, hit given on reverse strand
    • 20 = 4+16, read is unmapped, hit given on reverse strand but the read falls off the end of the chromosome.
  • Piping script:
for f in *.fastq
do
       /usr/local/bwa-0.5.9/bwa aln -l 28 -k 2 -t 16 -M 15 -I /media/disk-1/Dinh_Scratch/mm9Annotations-iGenomeUCSC/BWAIndex/mm9 $f > $f.sai
       /usr/local/bwa-0.5.9/bwa samse -n 1000 /media/disk-1/Dinh_Scratch/mm9Annotations-iGenomeUCSC/BWAIndex/mm9 $f.sai $f | ./parseIPtags-bwaSAM.pl $f.sam 
       
done;
for f in *.sam
do
       /usr/local/samtools-0.1.17/samtools view -uSt /media/disk-1/Dinh_Scratch/mm9Annotations-iGenomeUCSC/BWAIndex/genome.fai $f | /usr/local/samtools-0.1.17/samtools sort - sorted.$f
       /usr/local/samtools-0.1.17/samtools rmdup -S sorted.$f.bam rmdup.sorted.$f.bam
       /usr/local/samtools-0.1.17/samtools flagstat rmdup.sorted.$f.bam > rmdup.$f.LOG
       /usr/local/samtools-0.1.17/samtools flagstat sorted.$f.bam > sorted.$f.LOG
done;
File:ParseIPtags-bwaSAM.txt
  • RepeatTags (after rmdup)
5caC 171,839 (198,681)
5fC  177,485 (73,294)
IgG   12,904 (12,113)
Input 43,686 (36,396)
  • My original repeatTags files were really large, I have to figure out why samtools sort removed many lines of repeatTags:
 The number of lines in each file ( I checked that there are no empty lines printout by the perl script.
 19139320 YZ_UNC-5caC-DIP.fastq.repeatTags.sam
  4854616 YZ_UNC-5fC-DIP.fastq.repeatTags.sam
  4079042 YZ_UNC-IgG-DIP.fastq.repeatTags.sam
  6198485 YZ_UNC-Input.fastq.repeatTags.sam
  • UniqueTags (after rmdup)
5caC  22,864,806 (20,399,989)
5fC    4,147,394 (3,546,988)
IgG    8,010,605 (7,419,742)
Input 16,825,683 (16,295,873)

Balance unique hits[edit]

  • Use sort -R to get random hits
2611  samtools view rmdup.sorted.5fC.uniqueTags.bam > rmdup.sorted.5fC.uniqueTags.sam &
2615  samtools view rmdup.sorted.5caC.uniqueTags.bam | sort -R | head -3546988 > rmdup.sorted.5caC.uniqueTags.sam &
2616  samtools view rmdup.sorted.IgG.uniqueTags.bam | sort -R | head -3546988 > rmdup.sorted.IgG.uniqueTags.sam &
2617  samtools view rmdup.sorted.Input.uniqueTags.bam | sort -R | head -3546988 > rmdup.sorted.Input.uniqueTags.sam &

Peak calling of unique hits using MACS[edit]

macs14 -t rmdup.sorted.5caC.uniqueTags.sam -c rmdup.sorted.IgG.uniqueTags.sam -f SAM -g mm -n 5caC-vsIgG -w --call-subpeaks > macs14_5caCvsIgG_balanced &
macs14 -t rmdup.sorted.5caC.uniqueTags.sam -c rmdup.sorted.Input.uniqueTags.sam -f SAM -g mm -n 5caC-vsInput -w --call-subpeaks > macs14_5caCvsInput_balanced &
macs14 -t rmdup.sorted.5fC.uniqueTags.sam -c rmdup.sorted.IgG.uniqueTags.sam -f SAM -g mm -n 5fC-vsIgG -w --call-subpeaks > macs14_5fCvsIgG_balanced &
macs14 -t rmdup.sorted.5fC.uniqueTags.sam -c rmdup.sorted.Input.uniqueTags.sam -f SAM -g mm -n 5fC-vsInput -w --call-subpeaks > macs14_5fCvsInput_balanced &