Ns126:Encode Methylation

From ZhangLabWiki
Jump to navigation Jump to search

RRBS Data Analysis to Encode Project (Fastq)[edit]

RRBS Data Analysis to Encode Project (Bed/MethylFreq)[edit]

Data Download[edit]

Fastq and bed files can be downloaded from Encode Project. 101 RRBS data (bed files) were downloaded. 2,646,999 CpG loci were covered by 101 RRBS data while 866,979 CpG loci (32.8%) were detected in at least 80% samples.

  • File Download
cd /home/shg047/oasis/monod/rrbs_encode
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibMethylRrbs/files.txt
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibMethylRrbs/md5sum.txt
perl fastqDownloadRRBSEncode.pl files.txt
  • bed11 to bedGraph
cd /home/shg047/oasis/monod/rrbs_encode
for i in `ls *bed.gz`
do
zcat $i | grep -v "^track" |sort -k1,1 -k2,2n | awk '$5>9 {print $1"\t"$2"\t"$3"\t"$11}'> $i.bedGraph
done
  • Pearson Correlation
cd /home/shg047/oasis/monod/rrbs_encode

Aim 2: Haib39bioChain[edit]

  1. methylation haplotype block (MHB) calling with RRBS dataset

Background[edit]

  1. Basic:Single-end 40bp reads
  2. Encode: http://genome.ucsc.edu/ENCODE/downloads.html
  3. Encode|RRBS1: http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibMethylRrbs/
  4. Encode|RRBS2: http://genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=437674359_aUhx08DjchWwtBjyCYv61EB7Yy8S&c=chr1&g=wgEncodeHaibMethylRrbs
  5. Encode|Methy450K: http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibMethyl450/
  6. Encode|software: https://www.encodeproject.org/software/

Method[edit]

  1. Never Download Encode Data from UCSC. Encode Mainpage is great: https://www.encodeproject.org/search/?type=Experiment
  2. Download 39 biochain RRBS dataset from Encode Project (Fastq). Alignment with Bismark and merge all the BAM files
  3. make the haplotype calling with previous perl script. Finally MHB calling were down with Dr. Zhang perl script.
  4. Methyfreq based MHB calling were conducted with MethBed files download from UCSC.

Result[edit]

Summary Excel Haib Dataset manifest

  • Alignment: Maybe walltime time is too short , not all the samples were aligned completely, so I extend the walltime to 72 hours.
#PBS -q glean
#PBS -l nodes=1:ppn=8
#PBS -l walltime=72:00:00
bismark --bowtie2 --phred64-quals --fastq -L 30 -N 1 --multicore 2 /home/shg047/db/hg19/meth/bismark ../fastq_trim/HOT197_trimmed.fq.gz -o ../bam2
  • With above setting, the alignment time usage is about 16 hours
  • The time is propotional to the size of the fastq file
Sample N(reads) N(mapped) P(mapping) N(C) N(MCPG) N(MCHG) N(MCHH) N(UCPG) N(UCHG) N(UCHH) P(MCPG) P(MCHG) P(MCHH)
ENCFF000LVN 42665329 31372051 73.50% 227374577 9979920 566174 2568709 16516656 50531408 147211710 37.70% 1.10% 1.70%
ENCFF000LWL 37938401 25076517 66.10% 232580318 11488185 403775 1382726 44058434 51038500 124208698 20.70% 0.80% 1.10%
ENCFF000LVR 36779630 25226017 68.60% 216621540 12232832 438455 1846136 33041219 46480917 122581981 27.00% 0.90% 1.50%
ENCFF000LVW 36169867 22446797 62.10% 224690704 12356273 372366 871608 41484525 51598506 118007426 22.90% 0.70% 0.70%
ENCFF000LVB 35375019 22874896 64.70% 211937214 13517127 359167 1121727 31366464 47860817 117711912 30.10% 0.70% 0.90%
ENCFF000LUP 34214415 20212548 59.10% 221862905 16119238 340947 580628 49637185 51949984 103234923 24.50% 0.70% 0.60%
ENCFF000LWA 33924878 20895342 61.60% 226298397 14274674 319602 583229 52752860 53011846 105356186 21.30% 0.60% 0.60%
ENCFF000LUV 33798173 23173663 68.60% 249630261 10489904 225941 428583 62372541 57063134 119050158 14.40% 0.40% 0.40%
ENCFF000LWY 33715116 19999502 59.30% 213287381 13105506 316152 583822 46433586 48748868 104099447 22.00% 0.60% 0.60%
ENCFF000LWW 33642793 20282649 60.30% 216267059 15823449 340126 638093 43292443 49869268 106303680 26.80% 0.70% 0.60%
ENCFF000LWP 32858930 20654452 62.90% 219544115 13223974 319601 599593 48455363 50129553 106816031 21.40% 0.60% 0.60%
ENCFF000LUU 32815340 17866807 54.40% 196534662 12582324 421909 799252 48097172 45287260 89346745 20.70% 0.90% 0.90%
ENCFF000LVF 31276946 22130843 70.80% 241339808 11057578 348695 664774 64178792 55276052 109813917 14.70% 0.60% 0.60%
ENCFF000LXB 29681695 14796658 49.90% 160388755 8182288 279840 480680 41788931 37024968 72632048 16.40% 0.80% 0.70%
ENCFF000LWK 29349446 12860133 43.80% 137008343 8665157 302283 567461 29580229 31441673 66451540 22.70% 1.00% 0.80%
ENCFF000LWE 29211774 11162094 38.20% 115410335 8071486 277005 511615 23863725 26859483 55827021 25.30% 1.00% 0.90%
ENCFF000LVK 27335110 16938233 62.00% 193589572 12327753 322954 495345 51995628 43692466 84755426 19.20% 0.70% 0.60%
ENCFF000LVO 26190075 19276422 73.60% 150096025 7012551 300245 1332640 19959469 31819734 89671386 26.00% 0.90% 1.50%
ENCFF000LVA 25660341 16009728 62.40% 152217865 12748317 546338 1806399 26835875 36223859 74057077 32.20% 1.50% 2.40%
ENCFF000LWD 25467744 15459206 60.70% 149205383 11335949 535669 1829440 30226946 34383379 70894000 27.30% 1.50% 2.50%
ENCFF000LVU 23511285 15078397 64.10% 141659239 8631088 272002 734661 21915485 33328412 76777591 28.30% 0.80% 0.90%
ENCFF000LWO 22844877 14206467 62.20% 153005363 9813445 247506 461786 33745829 35223721 73513076 22.50% 0.70% 0.60%
ENCFF000LVI 22656883 13335953 58.90% 143634615 13714133 389619 650547 32241351 34451974 62186991 29.80% 1.10% 1.00%
ENCFF000LUQ 22247066 13650824 61.40% 150112632 9777330 287617 504659 38167142 34898963 66476921 20.40% 0.80% 0.80%
ENCFF000LUT 22240097 14849726 66.80% 162779228 9543204 285378 485629 41824570 38682843 71957604 18.60% 0.70% 0.70%
ENCFF000LWH 21943620 12872848 58.70% 137662143 11585659 414164 781429 29526002 32473470 62881419 28.20% 1.30% 1.20%
ENCFF000LVE 21473578 14519161 67.60% 157450743 7444292 221900 417095 40943018 36372568 72051870 15.40% 0.60% 0.60%
ENCFF000LWX 21208447 13218301 62.30% 145215571 9668358 223882 356827 34200783 33840207 66925514 22.00% 0.70% 0.50%
ENCFF000LVV 21080863 13197949 62.60% 142640926 10137634 255425 440789 32877645 33833167 65096266 23.60% 0.70% 0.70%
ENCFF000MLE 20338545 12207095 60.00% 110465581 11201397 457519 1827464 19226614 25845153 51907434 36.80% 1.70% 3.40%
ENCFF000LWS 20241908 11560661 57.10% 126245157 11387405 245297 424016 29669508 30429870 54089061 27.70% 0.80% 0.80%
ENCFF000LVZ 20058311 11004008 54.90% 120506134 9791631 242794 448332 25470114 28547035 56006228 27.80% 0.80% 0.80%
ENCFF000LWT 19407909 10887950 56.10% 119134111 9343415 203527 331921 25331949 28315688 55607611 26.90% 0.70% 0.60%
ENCFF000MLD 19184685 14418868 75.20% 162604971 6716952 305536 587608 54429999 38467374 62097502 11.00% 0.80% 0.90%
ENCFF000MLP 18102015 10032177 55.40% 109233809 10521402 215592 370895 25497472 26468867 46159581 29.20% 0.80% 0.80%
ENCFF000LVJ 17856591 11324104 63.40% 102317536 8591993 307833 1097735 17653908 24511494 50154573 32.70% 1.20% 2.10%
ENCFF000MLJ 16177211 8675428 53.60% 94678187 9430914 200777 318563 22135878 23265164 39326891 29.90% 0.90% 0.80%
ENCFF000LUN 10913594 3643074 33.40% 40450477 3011262 76408 126647 9141529 9549119 18545512 24.80% 0.80% 0.70%
ENCFF000MLM 5127779 2931330 57.20% 36675584 1281379 40172 75178 10725366 8687156 15866333 10.70% 0.50% 0.50%
  • MHB regions identification
#!/bin/csh
#PBS -N bam2MHB
#PBS -q pdafm
#PBS -l nodes=1:ppn=16
#PBS -l walltime=72:00:00
#PBS -o bam2MHB.log
#PBS -e bam2MHB.err
#PBS -V
#PBS -M shihcheng.guo@gmail.com
#PBS -m abe
#PBS -A k4zhang-group
cd /home/shg047/oasis/Haib/sortBam
# samtools cat -h header.sam -o haib.merge.bam *sort.bam
samtools sort -@ 16 haib.encode.merge.bam -o haib.merge.sort.bam
samtools index haib.merge.sort.bam
bedtools genomecov -bg -split -ibam haib.merge.sort.bam >   haib.merge.bam.pool.bed
awk '$4>9 { print $1"\t"$2"\t"$3}'  haib.merge.bam.pool.bed | bedtools merge -d 10 -i - > haib.RD10.genomecov.bed
awk '$3-$2>80 {print $1"\t"$2"\t"$3"\t"$3-$2+1}' haib.RD10.genomecov.bed > haib.RD10_80up.genomecov.bed
  • Statistic: haib.RD10_80up.genomecov.bed
cat haib.RD10_80up.genomecov.bed|awk '{sum+=$4} END { print "N = ", NR, "Sum = ", sum, " Average = ",sum/NR}' 
N =  120994 Sum =  17702479  Average =  146.309
  • haploinfo to MHB
cd /home/shg047/oasis/Haib/sortBam
/home/shg047/oasis/Haib/mhb/haib.RD10_80up.genomecov.bed
/home/shg047/oasis/Haib/mhb/haib.merge.sort.bam
/home/shg047/oasis/Haib/hapInfo2mld_blocks.pl
../mergedBam2hapInfo.pl ./haib.RD10_80up.genomecov.bed haib.merge.sort.bam > Haib.merge.RD10_80up.hapinfo.txt  # get hapinfo
../hapInfo2mld_block.pl ./Haib.merge.RD10_80up.hapinfo.txt 0.5 >  Haib.merge_RD10_80up.mld_blocks_r2-0.5.bed
  • MHB identified with different threshold: R-square from 0.1-0.9
/home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.1.bed
/home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.2.bed
/home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.3.bed
/home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.4.bed
/home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.5.bed
/home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.6.bed
/home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.7.bed
R-square threshold MHB counts
0.1 14933
0.2 13367
0.3 11667
0.4 9754
0.5 8155
0.6 7683
0.7 7445
bedtools intersect -wa -u /home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.5.bed
  • MHB calling based on RRBS Haib biochain data
cd /home/shg047/oasis/Haib/sortBam
samtools cat -h header.sam -o haib.encode.merge.bam *sort.bam
samtools sort haib.encode.merge.bam -o haib.encode.merge.sort.bam
samtools index haib.encode.merge.sort.bam
bedtools genomecov -bg -split -ibam haib.encode.merge.sort.bam >   haib.encode.merge.bam.pool.bed
awk '$4>9 { print $1"\t"$2"\t"$3}'  haib.encode.merge.bam.pool.bed | bedtools merge -d 10 -i - > haib.encode.RD10.genomecov.bed
awk '$3-$2>80 {print $1"\t"$2"\t"$3"\t"$3-$2+1}' haib.encode.RD10.genomecov.bed > haib.encode.RD10_80up.genomecov.bed
  • haploinfo to MHB
../hapInfo2mld_block.pl ./Haib.merge.RD10_80up.hapinfo.txt 0.5 >  Haib.merge_RD10_80up.mld_blocks_r2-0.5.bed