close
今天是生信星球陪你的第864天
大神一句話，菜鳥跑半年。我不是大神，但我可以縮短你走彎路的半年~
就像歌兒唱的那樣，如果你不知道該往哪兒走，就留在這學點生信好不好~
這裡有豆豆和花花的學習歷程，從新手到進階，生信路上有你有我!
文章標題：Characterization of an endoplasmic reticulum stress‐related signature to evaluate immune features and predict prognosis in glioma期刊：J Cell Mol Med影響因子:5.3雖然是個小小的5分文章，但涉及到的分析非常豐富，圖表也很多樣，我把他的數據拿來做例子進行分析，開啟一段更新咯。流程圖
構建模型部分的主要結果：
圖A是四個數據集中，用兩種算法（log-rank test和單因素cox）p<0.05 選出來的基因。B和C是lasso回歸的兩張經典圖表。
圖D是結合逐步回歸法，選出的16個基因構成的多因素cox模型。可以看到C-index值是0.86。這篇文章的後續用了兩個數據做測試集，效果都很好，那是因為一開始的時候篩選構建lasso回歸模型所使用的基因就是從四個數據集、兩種算法取出來的交集，這些基因都是與生存關係非常密切的，所以用他們構建出來的模型，在這四個數據集裡表現都不會差的，至於這算不算投機取巧嘛，無法定論，自己衡量就好。這個系列將會有好幾篇，今天先更新第一篇，整理數據。中間涉及到的一些算法和具體代碼的解析，需要一定的R語言基礎。如果你沒有系統學習過，可能看起來會比較困難。我的本職工作是做零基礎生信培訓的講師，可以拉到文末，有三個循環開課的培訓班的宣傳鏈接，可以來找我系統學習，學好基礎再上實戰就會絲滑啦。
1.數據整理的說明1.1 下載的數據四個數據的表達矩陣、生存信息表格。
TCGA的LGG和GBM是分開的。後續用代碼合併到一起。
1.2 表達矩陣的要求(1)就要是取過log、沒有異常值的矩陣，標準化過的也可以，因為不做差異分析。
(2)如果是轉錄組數據，要log之後的標準化數據。（因為這裡只做生存分析，不做差異分析，所以沒有必要找count數據）
(3)行名都轉換成gene symbol
(4)只要tumor樣本，不要normal
1.3 臨床信息的要求(1)要有每個樣本對應的生存時間和結局事件，並用代碼調整順序，與表達矩陣的樣本順序相同，嚴格一一對應。
(2)為了後續代碼統一，把生存時間和結局事件列名統一起來，都叫「time」和「event」。
(3)至於時間的單位，反正是分別計算的，統不統一無所謂。
1.4 關於數據來源的說明數據來源於哪個數據庫，不是主要矛盾。重點是分清楚是芯片數據還是轉錄組數據，是芯片的話要能找到探針注釋。如果是轉錄組數據，則需要分清楚是否標準化，進行了哪種標準化。反正是分別分析，又不需要合併，所以直接使用即可。
2.TCGA的RNA-seq數據Table S1顯示樣本數量691，TCGA的GBM數據173樣本，LGG有529樣本。
2.1 TCGA表達矩陣exp_gbm=read.table("TCGA-GBM.htseq_fpkm.tsv.gz",header=T,row.names=1,check.names=F)exp_lgg=read.table("TCGA-LGG.htseq_fpkm.tsv.gz",header=T,row.names=1,check.names=F)exp1=as.matrix(cbind(exp_gbm,exp_lgg))#行名轉換為genesymbollibrary(stringr)library(AnnoProbe)rownames(exp1)=str_split(rownames(exp1),"\\.",simplify=T)[,1];head(rownames(exp1))##[1]"ENSG00000242268""ENSG00000270112""ENSG00000167578""ENSG00000273842"##[5]"ENSG00000078237""ENSG00000146083"re=annoGene(rownames(exp1),ID_type="ENSEMBL");head(re)##SYMBOLbiotypesENSEMBLchrstart##1DDX11L1transcribed_unprocessed_pseudogeneENSG00000223972chr111869##2WASH7Punprocessed_pseudogeneENSG00000227232chr114404##3MIR6859-1miRNAENSG00000278267chr117369##4MIR1302-2HGlncRNAENSG00000243485chr129554##6FAM138AlncRNAENSG00000237613chr134554##7OR4G4Punprocessed_pseudogeneENSG00000268020chr152473##end##114409##229570##317436##431109##636081##753312library(tinyarray)exp1=trans_array(exp1,ids=re,from="ENSEMBL",to="SYMBOL")exp1[1:4,1:4]##TCGA-06-0878-01ATCGA-26-5135-01ATCGA-06-5859-01ATCGA-06-2563-01A##DDX11L10.00000000.02217550.000000000.01797799##WASH7P0.66739431.74299400.653346411.57566444##MIR6859-11.20067082.18805851.512291742.17610561##MIR1302-2HG0.00000000.00000000.025208570.00000000dim(exp1)##[1]56514702#刪除正常樣本Group=make_tcga_group(exp1)table(Group)##Group##normaltumor##5697exp1=exp1[,Group=="tumor"]#過濾表達量0值太多的基因exp1=exp1[apply(exp1,1,function(x){sum(x>0)>200}),]2.2 TCGA生存信息surv_gbm=readr::read_tsv("TCGA-GBM.survival.tsv")surv_gbm$TYPE="GBM"surv_lgg=readr::read_tsv("TCGA-LGG.survival.tsv")surv_lgg$TYPE="LGG"surv1=rbind(surv_gbm,surv_lgg)head(surv1)###Atibble:6x5##sampleOS`_PATIENT`OS.timeTYPE##<chr><dbl><chr><dbl><chr>##1TCGA-12-0657-01A1TCGA-12-06573GBM##2TCGA-32-1977-01A0TCGA-32-19773GBM##3TCGA-19-1791-01A0TCGA-19-17914GBM##4TCGA-28-1757-01A0TCGA-28-17574GBM##5TCGA-19-2624-01A1TCGA-19-26245GBM##6TCGA-41-4097-01A1TCGA-41-40976GBMtable(colnames(exp1)%in%surv1$sample)####FALSETRUE##6691s=intersect(colnames(exp1),surv1$sample)exp1=exp1[,s]surv1=surv1[match(s,surv1$sample),]colnames(surv1)[c(2,4)]=c("event","time")exp1[1:4,1:4]##TCGA-06-0878-01ATCGA-26-5135-01ATCGA-06-5859-01ATCGA-06-2563-01A##WASH7P0.667394251.742994000.653346411.57566444##MIR6859-11.200670782.188058481.512291742.17610561##AL627309.10.000000000.020662750.013869840.00000000##CICP270.096919810.010135270.000000000.02449211head(surv1)###Atibble:6x5##sampleevent`_PATIENT`timeTYPE##<chr><dbl><chr><dbl><chr>##1TCGA-06-0878-01A0TCGA-06-0878218GBM##2TCGA-26-5135-01A1TCGA-26-5135270GBM##3TCGA-06-5859-01A0TCGA-06-5859139GBM##4TCGA-06-2563-01A0TCGA-06-2563932GBM##5TCGA-41-2571-01A1TCGA-41-257126GBM##6TCGA-28-5207-01A1TCGA-28-5207343GBMidentical(colnames(exp1),surv1$sample)##[1]TRUE共691個有生存信息的tumor樣本。
3.CGGA芯片數據整理exp2=read.table("CGGA/CGGA.mRNA_array_301_gene_level.20200506.txt",header=T,row.names=1)surv2=read.table("CGGA/CGGA.mRNA_array_301_clinical.20200506.txt",sep="\t",header=T,check.names=F)head(surv2)##CGGA_IDTCGA_subtypesPRS_typeHistologyGradeGenderAgeOS##1CGGA_11ClassicalPrimaryGBMWHOIVFemale57155##2CGGA_124MesenchymalPrimaryGBMWHOIVMale53414##3CGGA_126MesenchymalPrimaryGBMWHOIVFemale501177##4CGGA_168MesenchymalPrimaryGBMWHOIVMale173086##5CGGA_172MesenchymalPrimaryGBMWHOIVFemale57462##6CGGA_195MesenchymalPrimaryGBMWHOIVMale48486##Censor(alive=0;dead=1)Radio_status(treated=1;un-treated=0)##111##211##311##401##511##611##Chemo_status(TMZtreated=1;un-treated=0)IDH_mutation_status##10Wildtype##20Wildtype##31Wildtype##41Wildtype##50Wildtype##61Wildtype##1p19q_Codeletion_statusMGMTp_methylation_status##1<NA>un-methylated##2<NA>un-methylated##3<NA>un-methylated##4<NA>un-methylated##5<NA>un-methylated##6<NA>un-methylateds=intersect(surv2$CGGA_ID,colnames(exp2))exp2=as.matrix(exp2[,s])surv2=surv2[match(s,surv2$CGGA_ID),]colnames(surv2)[c(9,8)]=c("event","time")exp2[1:4,1:4]##CGGA_11CGGA_124CGGA_126CGGA_168##A1BG-0.36036350.56495190.30477190.1749144##A1CF2.34136001.17079352.20815130.9652286##A2BP10.03451941.09640740.30248830.0949111##A2LD1-0.91305640.51088000.4253669-0.1949377head(surv2)##CGGA_IDTCGA_subtypesPRS_typeHistologyGradeGenderAgetimeevent##1CGGA_11ClassicalPrimaryGBMWHOIVFemale571551##2CGGA_124MesenchymalPrimaryGBMWHOIVMale534141##3CGGA_126MesenchymalPrimaryGBMWHOIVFemale5011771##4CGGA_168MesenchymalPrimaryGBMWHOIVMale1730860##5CGGA_172MesenchymalPrimaryGBMWHOIVFemale574621##6CGGA_195MesenchymalPrimaryGBMWHOIVMale484861##Radio_status(treated=1;un-treated=0)##11##21##31##41##51##61##Chemo_status(TMZtreated=1;un-treated=0)IDH_mutation_status##10Wildtype##20Wildtype##31Wildtype##41Wildtype##50Wildtype##61Wildtype##1p19q_Codeletion_statusMGMTp_methylation_status##1<NA>un-methylated##2<NA>un-methylated##3<NA>un-methylated##4<NA>un-methylated##5<NA>un-methylated##6<NA>un-methylatedidentical(surv2$CGGA_ID,colnames(exp2))##[1]TRUE我沒仔細看這個表達矩陣具體是哪種標準化，反正這個數據只是作為例子，我們就不深究到底是啥了，直接拿來用。注意如果是自己要用來發文章的數據，是要看清楚的，沒標準化需要自行標準化，cpm，tmp，fpkm都行。
5.GSE16011的整理library(tinyarray)geo=geo_download("GSE16011")library(stringr)#只要tumor樣本k=str_detect(geo$pd$title,"glioma");table(k)##k##FALSETRUE##8276geo$exp=geo$exp[,k]geo$pd=geo$pd[k,]#把行名轉換為genesymbolgpl=GEOquery::getGEO(filename="GPL8542_family.soft.gz",destdir=".")ids=gpl@dataTable@table[,1:2]library(clusterProfiler)library(org.Hs.eg.db)e2s=bitr(ids$ORF,fromType="ENTREZID",toType="SYMBOL",OrgDb="org.Hs.eg.db")ids=merge(ids,e2s,by.x="ORF",by.y="ENTREZID")ids=ids[,2:3]colnames(ids)=c("probe_id","symbol")exp4=trans_array(geo$exp,ids)surv4=geo$pd[,c(1,4,7,9)]exp4[1:4,1:4]##GSM405201GSM405202GSM405203GSM405204##A1BG7.2194016.0312616.7076816.330666##A2M12.20662613.26058712.35269313.044375##NAT15.1829897.2132726.5707336.343490##NAT25.2398425.2603735.6974354.859682head(surv4)##titleageatdiagnosis:ch1histology:ch1gender:ch1##GSM405201glioma844.57OD(gradeIII)Female##GSM405202glioma928.69OD(gradeIII)Male##GSM405203glioma1138.58OD(gradeIII)Male##GSM405204glioma1333.89OD(gradeIII)Male##GSM405205glioma2048.03OD(gradeIII)Male##GSM405206glioma2131.53OD(gradeIII)Maleidentical(rownames(surv4),colnames(exp4))##[1]TRUE臨床信息表格里沒有生存信息，所以從文章附件里找
從文章附件里提取臨床信息表格https://aacrjournals.org/cancerres/article/69/23/9065/553005/Intrinsic-Gene-Expression-Profiles-of-Gliomas-Are
包非常難裝，需要配置java，但這個需求不常用，直接跳過這個過程。
if(!file.exists("exp_surv4.Rdata")){library(tabulizer)f<-"00085472can092307-sup-stabs_1-6.pdf" re <- extract_tables(f,pages = 1:10)#提取前10頁的表格。str(re)re=do.call(rbind,re)re[1:4,1:4]colnames(re)=re[1,]re=re[-1,]re=data.frame(re)re[re==""]=NAlibrary(readr)re$Survival..years.=parse_double(re$Survival..years.,locale=locale(decimal_mark=","))re$Age.at.diagnosis=parse_double(re$Age.at.diagnosis,locale=locale(decimal_mark=","))dim(exp4)k=re$Reviewed.histological.diagnosis!="control";table(k)re=re[k,]re$Database.number=paste("glioma",re$Database.number)surv4$ID=rownames(surv4)surv4=merge(surv4,re,by.x="title",by.y="Database.number")colnames(surv4)[13:14]=c("event","time")table(surv4$event)surv4$time=as.integer(surv4$time*365)#天為單位surv4$event[surv4$event=="Lostto\rfollowup"]=NAtable(surv4$event)surv4$event=ifelse(surv4$event=="Alive",0,1)head(surv4)s=intersect(surv4$ID,colnames(exp4))exp4=exp4[,s]surv4=surv4[match(s,surv4$ID),]save(exp4,surv4,file="exp_surv4.Rdata")}load("exp_surv4.Rdata")exp4[1:4,1:4]##GSM405234GSM405235GSM405236GSM405237##A1BG7.6131506.9956667.3373936.287498##A2M13.16645713.00322513.81494212.991301##NAT15.5816727.4708947.8452106.088651##NAT25.2064644.6182454.8816525.011908head(surv4)##titleageatdiagnosis:ch1histology:ch1gender:ch1IDGender##1glioma10157.68GBM(gradeIV)FemaleGSM405234Female##2glioma10471.09GBM(gradeIV)MaleGSM405235Male##3glioma10552.20GBM(gradeIV)MaleGSM405236Male##4glioma10751.17GBM(gradeIV)MaleGSM405237Male##5glioma1138.58OD(gradeIII)MaleGSM405203Male##6glioma11137.25GBM(gradeIV)MaleGSM405238Male##Reviewed.histological.diagnosisAge.at.diagnosisKPSType.of.surgery##1GBM(gradeIV)57.6870Partialresection##2GBM(gradeIV)71.0980Partialresection##3GBM(gradeIV)52.2080Completeresection##4GBM(gradeIV)51.1770Stereotacticbiopsy##5OD(gradeIII)38.58<NA>Completeresection##6GBM(gradeIV)37.2580Stereotacticbiopsy##RadiotherapyChemotherapyeventtime##1yesno1226##2yesno176##3yesno1375##4yesno1149##5yesno13255##6yesno1343identical(surv4$ID,colnames(exp4))##[1]TRUE檢查和保存數據par(mfrow=c(2,2))boxplot(exp1[,1:10])boxplot(exp2[,1:10])boxplot(exp3[,1:10])boxplot(exp4[,1:10])
save(exp1,surv1,file="exp_surv1.Rdata")save(exp2,surv2,file="exp_surv2.Rdata")save(exp3,surv3,file="exp_surv3.Rdata")save(exp4,surv4,file="exp_surv4.Rdata")
全站熱搜
留言列表