Ns126:GMPOR: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Shicheng
>Shicheng
Line 219: Line 219:
  perl ~/bin/hapinfo2mhb.pl merge.txt.SumUniq 0.5 > mm9.MHB.0.5.txt  
  perl ~/bin/hapinfo2mhb.pl merge.txt.SumUniq 0.5 > mm9.MHB.0.5.txt  
  mkdir ../mhb
  mkdir ../mhb
  mv mm9.MHB.0.3.txt ../mhb
  mv mm9.MHB.0.5.txt ../mhb
  mv mm9.MHB.0.3.txt ../mhb
  mv mm9.MHB.0.3.txt ../mhb



Revision as of 01:39, 20 January 2017

Genome-wide DNA methylation profile of reprogramming

Method

  • Download human and mice reprogramming dataset
  • Merge genome-wide DNA methylation data of human and mice
  • Dendrogram based on genome-wide DNA methylation data(Figure 1a and 1b)
  • PCA analysis based on genome-wide DNA methylation data(Figure 2a and 2b)
  • heatmap based on DMR and DMRs in CpG Shore region (Figuare 1c,1d,1e and 1f)


wget ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1027nnn/GSM1027571/suppl/GSM1027571_DNA_CpG_coverage_E14_serum_LIF.bedGraph.gz &
wget ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1027nnn/GSM1027571/suppl/GSM1027571_DNA_CpG_methcounts_E14_serum_LIF.bedGraph.gz &
wget ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1027nnn/GSM1027572/suppl/GSM1027572_DNA_CpG_coverage_E14_2i_adapted.bedGraph.gz & 
wget ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1027nnn/GSM1027572/suppl/GSM1027572_DNA_CpG_methcounts_E14_2i_adapted.bedGraph.gz &
wget ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1127nnn/GSM1127953/suppl/GSM1127953_DNA_CpG_E14_serum_rep.bedGraph.gz &
wget ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1127nnn/GSM1127954/suppl/GSM1127954_DNA_CpG_E14_2i_adapted_rep.bedGraph.gz &

Result

  • MHB and entropy in scRNA-seq


Dendrogram and PCA analysis

  • PCA analysis based on genome-wide DNA methylation for human can be found here File:Figure2.PC13.pdf
  • PCA analysis based on genome-wide DNA methylation for mice can be found here

Differential Methylation Site Analysis

Distance between ES, IPS and HF

The distance between HF and ES as well as IPS is very large while actually IPS and ES are very similar, enclose figure showed ES, HF and two kinds of IPS (with two different inducing factor/method) in PC1-PC2 space.File:Ips ES HF.pdf.pdf, please check the code: File:PCAanalysis4GSE54115.R.txt

Relationship between Methylation linkage and Variation of Gene expression between ips and SCNT in mouse (mm9)

  • We don't have human data, so such analysis cannot be conducted in human
  • build analysis function so that Alice can use it. input: R2 and expression matrix.
  • input:
cd /media/Home_Raid1/zhl002/WGBS_mouse/
/media/Home_Raid1/zhl002/WGBS_mouse/iPSC_SCNT.R2Pvalue.high.txt

Mouse MHB

2016-9-7: Download 4 differentiation SRR files(SRR2011294/SRR2011297/SRR2011296/SRR2011295)

cd /home/shg047/oasis/mouse/fastq
perl ~/bin/trimGolrePBS.pl FastqMatchConfig.txt nonsubmit
qsub SRR2011294_1.fastq.gz.job
qsub SRR2011295_1.fastq.gz.job
qsub SRR2011296_1.fastq.gz.job
qsub SRR2011297_1.fastq.gz.job
perl ~/bin/bismark2bamPBS.pl FastqMatchConfig.txt nosubmit
qsub SRR2011294_1.fastq.gz.job
qsub SRR2011295_1.fastq.gz.job
qsub SRR2011296_1.fastq.gz.job
qsub SRR2011297_1.fastq.gz.job
  • RD90UP80
# bam merge
samtools merge -o output.bam $1 $2 $3
# bam split by CHR
for i in {1..19} X Y M
do
samtools view -b -q 20 Mouse.MHB.Alice.MergeBam.sort.bam chr$i > Mouse.MHB.Alice.chr$i.bam &
done
# bam to coverage
for i in `ls *bam`
do
bedtools genomecov -ibam $i -bg > merge.$i.bedcov
done
# Obtain D>=30 and Len>80 regions
for i in {1..19} X Y M
do
cd /home/shg047/oasis/mouse/RD
awk '$4>80{ print $1"\t"$2"\t"$3}' merge.chr$i.bam.gencov.bed | bedtools merge -d 80 -i - > merge.chr$i.bamRD90.bed
awk '$3-$2>90 {print $1"\t"$2"\t"$3"\t"$3-$2+1}' merge.chr$i.bamRD90.bed > merge.chr$i.bam.RD90_80up.bed
done
  • bam2hapinfo
perl /home/shg047/bin/bam2hapInfo2PBS.pl bam2hapinfo.config  submit bismark /home/shg047/oasis/db/mm9/mm9.chrom.sizes /home/shg047/oasis/db/mm9/MM9.CpG.positions.txt
cd /home/shg047/oasis/mouse/sortBam
perl ~/bin/bam2hapInfo2PBS.pl SamConfig.txt submit bismark /home/shg047/oasis/db/mm9/mm9.chrom.sizes /home/shg047/oasis/db/mm9/MM9.CpG.positions.txt
cd /home/shg047/oasis/mouse/alice/bam
perl ~/bin/bam2hapInfo2PBS.pl SamConfig.txt submit bisreadMapper /home/shg047/oasis/db/mm9/mm9.chrom.sizes /home/shg047/oasis/db/mm9/MM9.CpG.positions.txt
cd ../mergeBam/
perl ../saminfoPre4bam2hapinfo.pl > SamConfig.txt
perl /home/shg047/bin/bam2hapInfo2PBS.pl bam2hapinfo.config  submit bismark /home/shg047/oasis/db/mm9/mm9.chrom.sizes /home/shg047/oasis/db/mm9/MM9.CpG.positions.txt
qsub SRX1019864.job
qsub SRX1019865.job
qsub SRX1019866.job
qsub SRX1019867.job
  • merge haploinfo (cat) and Unique (hapinfoMerge.pl)
cat ************ | sort > Adult.sort.hapinfo.txt
cat ************ | sort > mESC.sort.hapinfo.txt
perl /home/shg047/oasis/mouse/mergeHapinfo/hapinfoMerge.pl Adult.sort.hapinfo.txt 
perl /home/shg047/oasis/mouse/mergeHapinfo/hapinfoMerge.pl mESC.sort.hapinfo.txt 
  • hapinfo to MHB
perl ~/bin/hapinfo2mhb.pl SRX202086.hapInfo.txt 0.1
  • Hapinfo2LD-R2(here: MHB regions,hapinfo2LDR2ByBed.pl)
cd /home/shg047/oasis/mouse/mergeHapinfo
perl ~/bin/hapinfo2BlocAvgR2.pl SRX1019864.hapInfo.txt > SRX1019864.R2.txt
cd /home/shg047/oasis/mouse/mergeHapinfo/mf
perl /home/shg047/bin/hapinfo2mf.pl /home/shg047/oasis/mouse/mergeHapinfo/mf > output.mf
File:4FD2.tm.png
  • hapinfo to R2 by chr10:126500749-126500761
perl ~/bin/hapinfo2LDR2.pl rlt chr10:126500749-126500761 < SRX202086.hapInfo.txt
  • hapinfo2mf (5me)
cd /home/shg047/oasis/mouse/mergeHapinfo/mf
perl ~/bin/hapinfo2mf.pl
grep chr10:100004267-100004288 output.mf
File:DD51.tm.png
  • Alice's own data
cd /home/shg047/oasis/mouse/alice/raw
ls *.txt.gz | awk -F. '{print $1}' | awk -F_ '{print $4}' | sort -u
  • Identify differential MHBs among ESC, Somatic tissue and IPS-SCNT groups
for i in {1..100}
do
perl ~/bin/randomSampleFromHaploInfo.pl tissue.hapinfo.txt > tissue.hapinfo.txt.$i
perl ~/bin/hapinfo2BlocAvgR2.pl tissue.hapinfo.txt.$i > R2.tissue.hapinfo.txt.$i
perl ~/bin/randomSampleFromHaploInfo.pl ipsNT.hapinfo.txt > ipsNT.hapinfo.txt.$i
perl ~/bin/hapinfo2BlocAvgR2.pl ipsNT.hapinfo.txt.$i > R2.ipsNT.hapinfo.txt.$i
perl ~/bin/randomSampleFromHaploInfo.pl ESC.hapinfo.txt > ESC.hapinfo.txt.$i
perl ~/bin/hapinfo2BlocAvgR2.pl ESC.hapinfo.txt.$i > R2.ESC.hapinfo.txt.$i
done


  • Differential Linkage MHBs Among IPS, STEM, TISSUE
data<-read.table("MatrixR2.txt",head=T)
x1<-grep("ESC",colnames(data))
x2<-grep("ips",colnames(data))
x3<-grep("tissue",colnames(data))
newdata<-data.frame(Y=c(rep("ESC",length(x1)),rep("IPS",length(x2)),rep("TISSUE",length(x3))),t(data),check.names=F)
Pvalue<-c()
for(i in 2:ncol(newdata)){
print(i)
Newdata<-na.omit(newdata[,c(1,i)])
if(length(unique(Newdata$Y))>1){
fit <- try(aov(Newdata[,2]~ Y, data=Newdata))
p<-summary(fit)1[1,5]
}else{
p<-NA  
}
Pvalue<-c(Pvalue,p)
}
rlt<-data.frame(MHB=colnames(newdata)[2:ncol(newdata)],Pvalue)
write.table(rlt,file="MHB.Pvalue.Bootstrap.txt",sep="\t",col.names=F,row.names=F,quote=F)
data1<-read.table("/oasis/tscc/scratch/shg047/mouse/hapinfo/group/MHB.Pvalue.Bootstrap.txt",head=F)
data2<-read.table("/home/shg047/oasis/mouse/hapinfo/tmp/MHB.R2.txt",head=T)
newdata<-data.frame(data2,Pvalue=data1[match(rownames(data2),data1[,1]),2])
write.table(newdata,file="MHB.R2.Pvalue.Bootstrap.txt",sep="\t",col.names=NA,row.names=T,quote=F)
data3<-read.table("/oasis/tscc/scratch/shg047/mouse/hapinfo/group/MHB.R2.Pvalue.Bootstrap.txt",head=F)
Subset1<-subset(data3,R2.ESC.hapinfo.txt<R2.ipsNT.hapinfo.txt & R2.tissue.hapinfo.txt<R2.ipsNT.hapinfo.txt)
dim(Subset1)
Subset2<-subset(data3,R2.ESC.hapinfo.txt<0.1 & R2.tissue.hapinfo.txt<0.1 & R2.ipsNT.hapinfo.txt>0.2)
dim(Subset2)
  • Two group wilcox test
data<-read.table("MatrixR2.txt",head=T)
x1<-grep("ESC",colnames(data))
x2<-grep("ips",colnames(data))
data<-data[,c(x1,x2)]
newdata<-data.frame(Y=c(rep("ESC",length(x1)),rep("IPS",length(x2))),t(data),check.names=F)
Pvalue<-c()
for(i in 2:ncol(newdata)){
 print(i)
 Newdata<-na.omit(newdata[,c(1,i)])
 if(length(unique(Newdata$Y))>1){
   fit <- try(wilcox.test(Newdata[,2]~ Y, data=Newdata))
   p<-fit$p.value
 }else{
   p<-NA  
 }
 Pvalue<-c(Pvalue,p)
}
rlt<-data.frame(MHB=colnames(newdata)[2:ncol(newdata)],Pvalue)
write.table(rlt,file="MHB.Pvalue.Bootstrap.txt",sep="\t",col.names=F,row.names=F,quote=F)
data1<-read.table("/oasis/tscc/scratch/shg047/mouse/hapinfo/group/MHB.Pvalue.Bootstrap.txt",head=F)
data2<-read.table("/home/shg047/oasis/mouse/hapinfo/tmp/MHB.R2.txt",head=T)
newdata<-data.frame(data2,Pvalue=data1[match(rownames(data2),data1[,1]),2])
write.table(newdata,file="MHB.R2.Pvalue.Bootstrap.txt",sep="\t",col.names=NA,row.names=T,quote=F)
data3<-read.table("/oasis/tscc/scratch/shg047/mouse/hapinfo/group/MHB.R2.Pvalue.Bootstrap.txt",head=F)
Subset1<-subset(data3,R2.ESC.hapinfo.txt<R2.ipsNT.hapinfo.txt & R2.tissue.hapinfo.txt<R2.ipsNT.hapinfo.txt)
dim(Subset1)
Subset2<-subset(data3,R2.ESC.hapinfo.txt<0.1 & R2.tissue.hapinfo.txt<0.1 & R2.ipsNT.hapinfo.txt>0.2)
dim(Subset2)

Human IPS/ESC MHB

  • Human hapinfo files, total 70. 1st part(59=49+10), 2nd part= 11 (ips and scnt)
cd /home/shg047/work/monod/hapinfo
ls *hapInfo.txt | grep -v CTT | grep -v PC | grep -v NC | grep -v LC | grep -v "\-7-" | grep -v 6- | grep -v "\-T-" | grep -v "\-P-" | grep -v merge | xargs -I {} cat {} >> merge.txt 
cd /home/shg047/work/Alice/human/mergeHapinfo
ls *hapinfo.txt | xargs -I {} cat {} >> /home/shg047/work/Alice/human/mergeHapinfo/merge.txt
# Go to GM since the memory of TSCC is not enough.
cd /home/shg047/work/monod/hapinfo
# already have the key, no need input the passwd
scp merge.txt shg047@genome-miner.ucsd.edu:/media/Home_Raid1/shg047/work/Alice/human/hapinfo
cd /media/Home_Raid1/shg047/work/Alice/human/hapinfo
perl ~/bin/hapinfoMerge.pl merge.txt 
perl ~/bin/hapinfo2mhb.pl merge.txt.SumUniq 0.3 > mm9.MHB.0.3.txt  
perl ~/bin/hapinfo2mhb.pl merge.txt.SumUniq 0.5 > mm9.MHB.0.5.txt 
mkdir ../mhb
mv mm9.MHB.0.5.txt ../mhb
mv mm9.MHB.0.3.txt ../mhb
/home/shg047/work/monod/hapinfo/MHB.0.5.txt
/home/shg047/work/monod/hapinfo/MHB.0.3.txt
  • 417 SRR download(complete, remember fastq-dump should be updated regular or else the download will be with lots of confusing mistake)
perl srrdownload.pl
  • 417 Alignment by bismark and BAM (complete)
cd /home/shg047/work/Alice/human
perl bismarkAutoSingle.pl submit
  • 417 sortBam(complete)
cd /home/shg047/work/Alice/bam
perl ~/bin/bam2sortBam.pl sortc submit
  • Remove cancer samples
  • 417 hapinfo (complete)
cd /home/shg047/work/Alice/sortBam
perl ~/bin/SaminfoPre4hapinfo.pl > sampleconfig.txt
perl ~/bin/bam2hapInfo2PBS.pl sampleconfig.txt non bismark /home/shg047/oasis/db/hg19/hg19.chrom.sizes /home/shg047/oasis/db/hg19/HsGenome19.CpG.positions.txt
  • coverage, R2, methylation level
WNT5A chr14:29318659-29338701
PPP3CA chr3:136333030-136600691
sort the new bam files

Conclusion

  • For human ips and SCNT data, the methylation status of 4377908 CpG were shared by 11 samples(2 ES, 1 Fibroblast, 4 ips and 4 SCNT). PCA analysis based on the genome-wide DNA methylation data showed the first and second principle component can explain almost 50% variances. In the dimension of first component fibroblast were obviously separated with ES,ips and SCNT. In the dimension of the second component, ES was separated with ips and SCNT with clear boundary.
  • For mice ips and SCNT data, the methylation status of 4909250 CpG were shared by 9 samples(1 ES, 4 ips and 4 SCNT). PCA analysis based on the genome-wide DNA methylation data showed the first and second principle component can explain almost 50% variances. In the dimension of first component fibroblast were obviously separated with ES,ips and SCNT. In the dimension of the second component, ES was separated with ips and SCNT with clear boundary.

Highlight

Human ES/iPS Dataset

GSM1385973 SCNT-1
GSM1385974 HDF
GSM1385975 HESO-7
GSM1385976 HESO-8
GSM1385977 SCNT-2
GSM1385978 SCNT-3
GSM1385979 SCNT-4
GSM1385980 iPS-S1
GSM1385981 iPS-S2
GSM1385982 iPS-R2
GSM1385983 iPS-R1
SCNT1 SCNT2 SCNT3 SCNT4 iPSC_S1 iPSC_S2 iPSC_R1 iPSC_R2
SRX542415 SRX542419 SRX542420 SRX542421 SRX542422 SRX542423 SRX542425 SRX542424
SRR1286295 SRR1286459 SRR1286502 SRR1286538 SRR1286579 SRR1286620 SRR1286703 SRR1286661
SRR1286296 SRR1286460 SRR1286503 SRR1286539 SRR1286580 SRR1286621 SRR1286704 SRR1286662
SRR1286297 SRR1286461 SRR1286504 SRR1286540 SRR1286581 SRR1286622 SRR1286705 SRR1286663
SRR1286298 SRR1286462 SRR1286505 SRR1286541 SRR1286582 SRR1286623 SRR1286706 SRR1286664
SRR1286299 SRR1286463 SRR1286506 SRR1286542 SRR1286583 SRR1286624 SRR1286707 SRR1286665
SRR1286300 SRR1286464 SRR1286507 SRR1286543 SRR1286584 SRR1286625 SRR1286708 SRR1286666
SRR1286301 SRR1286465 SRR1286508 SRR1286544 SRR1286585 SRR1286626 SRR1286709 SRR1286667
SRR1286302 SRR1286466 SRR1286509 SRR1286545 SRR1286586 SRR1286627 SRR1286710 SRR1286668
SRR1286303 SRR1286467 SRR1286510 SRR1286546 SRR1286587 SRR1286628 SRR1286711 SRR1286669
SRR1286304 SRR1286468 SRR1286511 SRR1286547 SRR1286588 SRR1286629 SRR1286712 SRR1286670
SRR1286305 SRR1286469 SRR1286512 SRR1286548 SRR1286589 SRR1286630 SRR1286713 SRR1286671
SRR1286306 SRR1286470 SRR1286513 SRR1286549 SRR1286590 SRR1286631 SRR1286714 SRR1286672
SRR1286307 SRR1286471 SRR1286514 SRR1286550 SRR1286591 SRR1286632 SRR1286715 SRR1286673
SRR1286308 SRR1286472 SRR1286515 SRR1286551 SRR1286592 SRR1286633 SRR1286716 SRR1286674
SRR1286309 SRR1286473 SRR1286516 SRR1286552 SRR1286593 SRR1286634 SRR1286717 SRR1286675
SRR1286310 SRR1286474 SRR1286517 SRR1286553 SRR1286594 SRR1286635 SRR1286718 SRR1286676
SRR1286311 SRR1286475 SRR1286518 SRR1286554 SRR1286595 SRR1286636 SRR1286719 SRR1286677
SRR1286312 SRR1286476 SRR1286519 SRR1286555 SRR1286596 SRR1286637 SRR1286720 SRR1286678
SRR1286313 SRR1286477 SRR1286520 SRR1286556 SRR1286597 SRR1286638 SRR1286721 SRR1286679
SRR1286314 SRR1286478 SRR1286521 SRR1286557 SRR1286598 SRR1286639 SRR1286722 SRR1286680
SRR1286315 SRR1286479 SRR1286522 SRR1286558 SRR1286599 SRR1286640 SRR1286723 SRR1286681
SRR1286316 SRR1286480 SRR1286523 SRR1286559 SRR1286600 SRR1286641 SRR1286724 SRR1286682
SRR1286317 SRR1286481 SRR1286524 SRR1286560 SRR1286601 SRR1286642 SRR1286725 SRR1286683
SRR1286318 SRR1286482 SRR1286525 SRR1286561 SRR1286602 SRR1286643 SRR1286726 SRR1286684
SRR1286319 SRR1286483 SRR1286526 SRR1286562 SRR1286603 SRR1286644 SRR1286727 SRR1286685
SRR1286320 SRR1286484 SRR1286527 SRR1286563 SRR1286604 SRR1286645 SRR1286728 SRR1286686
SRR1286321 SRR1286485 SRR1286528 SRR1286564 SRR1286605 SRR1286646 SRR1286729 SRR1286687
SRR1286322 SRR1286486 SRR1286529 SRR1286565 SRR1286606 SRR1286647 SRR1286730 SRR1286688
SRR1286323 SRR1286487 SRR1286530 SRR1286566 SRR1286607 SRR1286648 SRR1286731 SRR1286689
SRR1286324 SRR1286488 SRR1286531 SRR1286567 SRR1286608 SRR1286649 SRR1286732 SRR1286690
SRR1286325 SRR1286489 SRR1286532 SRR1286568 SRR1286609 SRR1286650 SRR1286733 SRR1286691
SRR1286326 SRR1286490 SRR1286533 SRR1286569 SRR1286610 SRR1286651 SRR1286734 SRR1286692
SRR1286327 SRR1286491 SRR1286534 SRR1286570 SRR1286611 SRR1286652 SRR1286735 SRR1286693
SRR1286328 SRR1286492 SRR1286535 SRR1286571 SRR1286612 SRR1286653 SRR1286736 SRR1286694
SRR1286329 SRR1286493 SRR1286536 SRR1286572 SRR1286613 SRR1286654 SRR1286737 SRR1286695
SRR1286330 SRR1286494 SRR1286537 SRR1286573 SRR1286614 SRR1286655 SRR1286738 SRR1286696
SRR1286331 SRR1286495 SRR1286574 SRR1286615 SRR1286656 SRR1286739 SRR1286697
SRR1286332 SRR1286496 SRR1286575 SRR1286616 SRR1286657 SRR1286740 SRR1286698
SRR1286333 SRR1286497 SRR1286576 SRR1286617 SRR1286658 SRR1286741 SRR1286699
SRR1286334 SRR1286498 SRR1286577 SRR1286618 SRR1286659 SRR1286742 SRR1286700
SRR1286335 SRR1286499 SRR1286578 SRR1286619 SRR1286660 SRR1286743 SRR1286701
SRR1286336 SRR1286500 SRR1286702
SRR1286337 SRR1286501

Supplemenatry

  • For mice

Reference

  • 1.Ma H, Morey R, O'Neil RC, He Y, Daughtry B, Schultz MD, et al. Abnormalities in human pluripotent cells due to reprogramming mechanisms. Nature. 2014;511(7508):177-83. Epub 2014/07/11.
  • 2.
  • 3.