Ns126:Calendar/NOTES/2015-4-21

From ZhangLabWiki
Revision as of 03:35, 22 April 2015 by >Shicheng
Jump to navigation Jump to search
  • Now, I do know how to insert figures, but I don't know how to insert tables into MediaWiki. Therefore, I write a small code to transfer routine table to wikitable.
  • I must say: mediawiki is really a excellent tool to manage the lab and monitor the progress or task to help us complete the project quickly.
perl table2wikitable.pl input.table.txt
#!/usr/bin/perl
#table2wikitable.pl
use strict;
use Cwd;
chdir getcwd;

my $input=@ARGV[0];

open F,$input;
print "{| class=\"wikitable\" style=\"text-align: right; color: red;\"\n";
while(<F>){
chomp;
my @line=split /\t/;
my $tmp=join("||",@line);
print "| $tmp\n|-\n";
}
print"|}\n";
 
  • corresponding table is as the following. Pan-Cancer Methylation 450K dataset were collected from TCGA Project to identify optional padloc regions.
Symbol Case Control Cancer 肿瘤名称
KIRC 160 160 Kidney renal clear cell carcinoma 肾透明细胞癌
BRCA 92 92 Breast invasive carcinoma 浸润性乳腺癌
THCA 56 56 Thyroid carcinoma 甲状腺癌
HNSC 50 50 Head and Neck squamous cell carcinoma 头颈部鳞状细胞癌
PRAD 49 49 Prostate adenocarcinoma 前列腺癌
LIHC 49 49 Liver hepatocellular carcinoma 肝癌
KIRP 45 45 Kidney renal papillary cell carcinoma 肾乳头状细胞癌
LUSC 41 41 Lung squamous cell carcinoma 肺腺癌
COAD 39 39 Colon adenocarcinoma 结肠癌
UCEC 30 30 Uterine Corpus Endometrial carcinoma 子宫内膜癌
LUAD 26 26 Lung adenocarcinoma 肺鳞癌
  • In order to publish MONOD paper to a good paper quickly. 1) fast 2) large data
    • collect methylation 450 data of cancer tissues and normals from TCGA Project
      • sample size
      • clinical information: gender, age
tar xvf 08177614-f305-4fbf-84ca-fd2fbfe26755.tar   #
tar xvf 548dd1cf-84a1-4d96-86f0-b47d6300daa6.tar   #
tar xvf 9cd95b1c-782c-478d-9ec9-de7c9c441cc4.tar   #
tar xvf c5da6cd3-0266-4ebb-bd11-713ffe9b5ef9.tar   #
  • focus on 3 cancers which we have preliminary data (lung, colon, pancreatic)
library("stringr")
for(cancer in c("COAD","LUAD","LUSC","PAAD")){
dir<-paste("/home/sguo/monod/data/",tolower(cancer),"/DNA_Methylation/JHU_USC__HumanMethylation450/Level_3",sep="")
setwd(dir)
pattern=paste("jhu-usc.edu_",cancer,".*",sep="")
print (pattern)
file=list.files(pattern=pattern)
idv<-unique(as.array(str_extract(file,"TCGA-[0-9|a-z|A-Z]*-[0-9|a-z|A-Z]*")))
pairidv<-c()
for (i in 1:length(idv)){
  t1<-paste(idv[i],"-01",sep="") 
  t2<-paste(idv[i],"-11",sep="")
  if(all(any(grepl(t1,file)),any(grepl(t2,file)))){
    pairidv<-c(pairidv,t1,t2)
  }
}
l1<-length(pairidv)
l2<-length(file)
id1<-lapply(lapply(strsplit(file,"[.]"),function(x) x[6]),function(x) substr(x,1,15))
id2<-lapply(lapply(strsplit(file,"[.]"),function(x) x[6]),function(x) substr(x,14,15))
sam<-lapply(lapply(strsplit(file,"[.]"),function(x) x[6]),function(x) substr(x,1,15))
tab<-table(unlist(lapply(lapply(strsplit(file,"[.]"),function(x) x[6]),function(x) substr(x,14,15))))
c1<-tab[which(names(tab)=="01")]
c2<-tab[which(names(tab)=="11")]
c3<-length(pairidv)
print(c((c1),(c2),Pair=c3))
}
  • The sample size for the genome-wide DNA methyaltion dataset for 3 cancers in TCGA Project are as the following:
Cancer Normal Paired Total
COAD 312 38 76 350
LUAD 473 32 58 505
LUSC 370 42 80 412
PAAD 184 10 20 194
Sum 1339 122 234 1461
    • totally, 1339 cancer tissues and 122 normal tissues were used to refine the padloc regions.
    • firstly, we should check the methylation status of the CpG sites/regions in the normal tissues.
1) hypermethylation in cancer/patient plasma
2) hypomethylation in normal/health plasma
3) large methylation difference (delta beta)
4) high group specificity index (GSI)
5) high ratio release to plasma

File:GSI2.png

  • s(j): average methylation in individual group j
  • S(max): average methylation in the group with highest methylation level
# group specificity index,GSI
setwd("../pan")
load("PanPairMethData.RData")
pheno=data$pheno
xmean <- rowsum(data[,2:ncol(data)], pheno)/table(pheno)
gsi<-apply(xmean,2,function(x) (length(x)-sum(x)/max(x))/(length(x)-1))
rlt<-data.frame(pheno=names(table(pheno)),xmean[,match(names(sort(gsi,decreasing=T)[1:5]),colnames(xmean))])
write.table(rlt,file="pancancer.high.gsi.site.txt",sep="\t",quote=F,col.names=NA,row.names=T)
pdf("gsi.distribution.pdf")
hist(gsi,main="Histogram of Group Specific Index")
dev.off()
num90<-sum(gsi>0.85)
xmean90<-xmean[,match(names(sort(gsi,decreasing=T)[1:num90]),colnames(xmean))]
xx1<-table(data$pheno)[as.numeric(names(table(apply(xmean90,2,function(x) which.max(x)))))]
xx2<-table(apply(xmean90,2,function(x) which.max(x)))
rlt2<-data.frame(sample_size=xx1,xx2)
write.table(rlt2,file="target.status.specifi.txt",sep="\t",quote=F,col.names=NA,row.names=T)
  • run the code and you can find the strongest group specific index CpGs with TCGA Pancancer dataset.
Sample_size cg16579555 cg26240185 cg12587766 cg15375239 cg10157975
BRCA-01 0.01 0.014 0.021 0.012 0.011
BRCA-11 0.01 0.011 0.02 0.011 0.011
COAD-01 0.011 0.012 0.605 0.011 0.353
COAD-11 0.011 0.015 0.033 0.012 0.022
HNSC-01 0.011 0.012 0.021 0.013 0.031
HNSC-11 0.011 0.011 0.022 0.012 0.012
KIRC-01 0.009 0.01 0.016 0.019 0.011
KIRC-11 0.009 0.009 0.016 0.01 0.009
KIRP-01 0.01 0.026 0.021 0.016 0.01
KIRP-11 0.01 0.011 0.022 0.012 0.011
LIHC-01 0.396 0.498 0.042 0.468 0.028
LIHC-11 0.057 0.099 0.024 0.133 0.015
LUAD-01 0.01 0.014 0.022 0.012 0.012
LUAD-11 0.011 0.013 0.024 0.013 0.013
LUSC-01 0.009 0.01 0.017 0.01 0.009
LUSC-11 0.008 0.01 0.018 0.011 0.01
PRAD-01 0.011 0.011 0.02 0.012 0.011
PRAD-11 0.01 0.012 0.019 0.012 0.011
THCA-01 0.011 0.013 0.028 0.013 0.013
THCA-11 0.011 0.013 0.024 0.012 0.013
UCEC-01 0.013 0.014 0.02 0.012 0.026
UCEC-11 0.011 0.015 0.022 0.012 0.014
  • clearly, we can find cancer specific hypermethylated genes with previous method. Then, how about the distribution of the GSI.

File:Gsi.distribution.jpeg

  • 745 CpGs or regions whose GSI>0.9. these CpGs/regions showed high tissue/status specific.
  • 3181 CpGs or regions whose GSI>0.85. these CpGs/regions showed high tissue/status specific.
  • The question is that whether they can be found in patients plasma? (Dr.Zhang might have some information)
  • 3181 tissue and status specific hypermethylation CpG sites were showed as the following table.
Sample_size Var1 Freq
BRCA-01 92 1 33
BRCA-11 92 2 1
COAD-01 39 3 1695
COAD-11 39 4 16
HNSC-01 50 5 143
KIRC-01 160 7 1
KIRC-11 160 8 1
KIRP-01 45 9 20
KIRP-11 45 10 1
LIHC-01 49 11 768
LIHC-11 49 12 40
LUAD-01 26 13 1
LUSC-01 41 15 8
PRAD-01 49 17 271
PRAD-11 49 18 1
THCA-01 56 19 4
THCA-11 56 20 7
UCEC-01 30 21 168
UCEC-11 30 22 2
  • Here, I choose top 5 regions to represent such specific tissue and it's specific status (cancer or normal), code is as the following:
gsi2<-gsi[order(gsi,decreasing=T)]
rlt3<-c()
for(i in seq(1,22,by=2)){
  z<-0
  j<-0
  while(j<length(gsi2) & z<5){
  j=j+1
  max<-which.max(xmean[,match(names(gsi2[j]),colnames(xmean))])
  if(max==i){
    z<-z+1
    rlt3<-rbind(rlt3,(c(names(gsi2[j]),names(table(pheno))[i])))
  }
}
}
rlt3<-data.frame(rlt3)
inf<-read.table("/home/sguo/monod/data/paad/DNA_Methylation/JHU_USC__HumanMethylation450/Level_3/jhu-usc.edu_PAAD.HumanMethylation450.9.lvl-3.TCGA-S4-A8RM-01A-11D-A378-05.txt",head=F,skip=2,sep="\t")
cpg<-rlt3[,1]
cancer<-substr(rlt3[,2],1,4)
gene<-inf[match(rlt3[,1],inf[,1]),3]
chr<-inf[match(rlt3[,1],inf[,1]),4]
start<-inf[match(rlt3[,1],inf[,1]),5]-60
end<-inf[match(rlt3[,1],inf[,1]),5]+60
GSI<-gsi[match(rlt3[,1],names(gsi))]
rlt4<-data.frame(cpg,cancer,GSI,gene,chr,start,end)
write.table(rlt4,file="target.status.specifi.site.txt",sep="\t",quote=F,col.names=NA,row.names=T)
  • The results were showed as the following table. Top 5 important biomarkers for each cancer of the specific tissue and status were list.
cpg cancer GSI gene chr start end
cg18565473 BRCA 0.932682778 ETS1 11 128392042 128392162
cg23884187 BRCA 0.9157933441 C20orf95 20 37275009 37275129
cg14052221 BRCA 0.912760996 PSAT1 9 80911998 80912118
cg24797187 BRCA 0.9120597976 AFF3 2 100175708 100175828
cg18943693 BRCA 0.9116756297 1 155043501 155043621
cg12587766 COAD 0.9628664076 LIFR 5 38556375 38556495
cg10157975 COAD 0.9590682334 ZNF304 19 57862382 57862502
cg23977631 COAD 0.9547387465 LONRF2 2 100938739 100938859
cg04117229 COAD 0.9529809636 SPG20 13 36920753 36920873
cg09854653 COAD 0.9527461786 QKI 6 163834843 163834963
cg03988778 HNSC 0.9412487691 SVIP 11 22850831 22850951
cg08211306 HNSC 0.9364418236 ENPP4 6 46097724 46097844
cg26968387 HNSC 0.9325520916 ZNF420 19 37569208 37569328
cg03280624 HNSC 0.9251408321 ZNF583 19 56915595 56915715
cg00471966 HNSC 0.9237859158 ZNF420 19 37569290 37569410
cg26228351 KIRC 0.8668214436 KIF21B 1 200992596 200992716
cg00593900 KIRC 0.8404290268 ANGPTL6 19 10206686 10206806
cg11697226 KIRC 0.8393271691 TNFRSF11A 18 59992325 59992445
cg09865339 KIRC 0.83639651 GPC2;STAG3 7 99774875 99774995
cg02632185 KIRC 0.8242621029 MAST4 5 66299726 66299846
cg15598442 KIRP 0.9265204439 1 25175003 25175123
cg26622232 KIRP 0.9074060913 OXR1 8 107669727 107669847
cg16326979 KIRP 0.8986303106 OXR1 8 107670101 107670221
cg17031478 KIRP 0.8970812881 HOXC4;HOXC5 12 54427113 54427233
cg17136799 KIRP 0.8884003413 OXR1 8 107669723 107669843
cg16579555 LIHC 0.9680462923 RNF135 17 29298292 29298412
cg26240185 LIHC 0.966346469 FAR1 11 13690097 13690217
cg15375239 LIHC 0.961490202 SPINT2 19 38755227 38755347
cg15969216 LIHC 0.9580311599 TSC22D1 13 45150202 45150322
cg03326059 LIHC 0.9580238498 FAR1 11 13690100 13690220
cg13215643 LUAD 0.8553268659 DACT1 14 59104765 59104885
cg07017994 LUAD 0.8425613476 EPHB6 7 142552854 142552974
cg12487147 LUAD 0.8391970226 HSD17B8 6 33172382 33172502
cg26615830 LUAD 0.8371456733 MSX1 4 4861270 4861390
cg21929771 LUAD 0.8207003816 PTPRU 1 29586520 29586640
cg07240673 LUSC 0.8990393087 CLUAP1 16 3550848 3550968
cg18772127 LUSC 0.8862042556 CMTM7 3 32443436 32443556
cg08562243 LUSC 0.8748813166 CLUAP1 16 3551109 3551229
cg02566698 LUSC 0.8726750003 CLUAP1 16 3550872 3550992
cg02379764 LUSC 0.8697737274 CLUAP1 16 3550968 3551088
cg07635623 PRAD 0.9556665683 SERPINB1 6 2841815 2841935
cg14283569 PRAD 0.9475569263 19 51416153 51416273
cg05098590 PRAD 0.9434078407 ADD3 10 111767319 111767439
cg10938374 PRAD 0.942601773 IER3 6 30711998 30712118
cg16232979 PRAD 0.9424310073 TPM4 19 16187571 16187691
cg27115721 THCA 0.8819506898 FAM49A 2 16790278 16790398
cg04358131 THCA 0.8775891556 MAFK 7 1572192 1572312
cg10540754 THCA 0.8609459396 FAM49A 2 16790310 16790430
cg22749810 THCA 0.8552029028 RNF213 17 78237329 78237449
cg12822074 THCA 0.8497825469 RTN4RL2 11 57243805 57243925
cg03221247 UCEC 0.9555805745 LYPLAL1 1 219347398 219347518
cg15494117 UCEC 0.9514218403 TERC 3 169482835 169482955
cg15599946 UCEC 0.9474897559 TERC 3 169482839 169482959
cg18985581 UCEC 0.946277345 14 105512153 105512273
cg02665570 UCEC 0.9458580151 LYPLAL1 1 219347280 219347400