已测试:TCGA差异表达分析和生存分析
TCGA数据库介绍请自行查阅资料,推荐链接:https://www.jianshu.com/p/dd066035b399

上图可以看到,TCGA的数据包括了基因测序,基因拷贝数分析,甲基化分析,转录组分析,miRNA分析,临床数据。
那么,我们可以分析哪些数据呢?
基因测序数据,我们可以看出基因突变;基因拷贝数,我们可以看到基因组中哪些基因拷贝数多;甲基化数据可以找出差异的甲基化位点或者区域;转录组可以知道哪些基因高表达,哪些低表达;miRNA分析,可以找出差异的miRNA,或者找出跟癌症相关的miRNA。它还有临床的数据,基因的数据可以和临床数据结合,比如不同的分组,是否突变不同,对生存期是否不同。想知道某个基因的表达,对病人的预后是否影响,可以把转录组和预后进行分析。分子水平也可以联合分析,如转录组可以和miRNA联合起来,可以看出哪些miRNA会和基因的表达有关,哪个miRNA会对某个基因有调控的作用。甲基化也可以和转录组联合,甲基化的水平对基因表达的影响。

上图是样品的名称,sample中,01-09就是肿瘤组织,10-19是正常组织,20-29是正常人。
进入TCGA官网
https://portal.gdc.cancer.gov/

左上角的选择case, 可以选择不同组织,根据左侧列表进行筛选。我们选择前列腺癌,如下图,右边显示的是2326个数据。

下面program里面选择TCGA。如果我们要选择多个数据库,也可以点选。
case里面选择完了,我们继续选择file里面的数据。

我们选择转录本的数据。

加入购物车以后,点击购物车,则可以看到如下可以下载的页面,我们可以点击各个下载的地方,如download里面,有一个manifest文件,是需要直接下载的,此文件是解释文件,解释我们即将下载的文件。cart是基因表达量的文件。clinical是临床数据。biospecimen人群数据。filemetadata是名单库,需要和我们的下载的名单对应起来。
以上文件依次下载,但是cart文件可能经常下载不了,网络问题等导致下载失败,我们就需要用到命令行下载。
需要先进行工具下载,TCGA主页右上角有一个GDC apps,下载data transfer tool。根据自己的系统下载,不需要安装,因为软件是二进制的。

找到manifest文件,放在同一个目录下,打开dos环境,运行gdc-client.exe download -m MANIFEST.txt
一定下载最新版本的GDC工具,不然会报错,TCGA经常更新。


上图显示正在下载。
下一步则进行数据整理,整理成矩阵以后,则进行id转为gene symbol,整理生存数据。
把所有的文件都放在一个目录下。这个可以手动,也可以通过代码。
用perl移动文件至同一个目录,这里要注意,不能把miRNA和mRNA的数据下载到同一个目录下,不然都给融合到一个文件夹了。具体的文件可下载附件,或者留言索取。如下图所示。

进入cmd,进入目录文件夹,输入以下命令,回车,就可以了。
C:\Users\Administrator\Desktop\111>perl putFilesToOneDir.pl
就会出现以下一个文件夹,就是完全整理到一个文件夹了。

然后全部解压。


然后把mRNA的数据整合到一起,用mRNAmerge脚本和metadata文件,前面已经下载过,复制到上面解压后文件的文件夹。同样按照perl的操作方式,进入到目录以后,运行perl,空格,后面跟着mRNAmerge脚本,再加上metadata文件。
C:\Users\Administrator\Desktop\111\files>perl mRNA_merge.pl metadata.cart.2020-02-20.json
回车后,这个运行需要一定的时间。

运行完了以后,系统会告诉你,normal有多少,tumor有多少。这个数字要记录好,后面会用得到。
运行完毕,会出现一个mRNAmatrix.txt文件,就是我们后续分析需要用到的。

打开这个文件,我们可以看到,前51个就是normal。关键的字段是11,而不是1到9。

下一步,我们要把ensemble ID转化为基因名。

对应文件需要到ensemble的官网去下载。http://asia.ensembl.org/index.html

download,下载。



下载好以后解压,然后运行perl ensemblToSymbol.pl Homo_sapiens.GRCh38.99.chr.gtf mRNAmatrix.txt mRNA.symbol.txt(这里是我们自己命名的输出文件名字。)
运行后,则基因的名字转换完毕。

foldChange=1 #设置差异倍数
padj=0.05 #设置P值
library("edgeR")
rt=read.table("mRNA.symbol.txt",sep="\t",header=T,check.names=F) #改成自己的文件名
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data=data[rowMeans(data)>1,]
group=c(rep("normal",51),rep("tumor",489)) #按照癌症和正常样品数目修改
design <- model.matrix(~group)
y <- DGEList(counts=data,group=group)
y <- calcNormFactors(y)
y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)
et <- exactTest(y,pair = c("normal","tumor"))
topTags(et)
ordered_tags <- topTags(et, n=100000)
allDiff=ordered_tags$table
allDiff=allDiff[is.na(allDiff$FDR)==FALSE,]
diff=allDiff
newData=y$pseudo.counts #差异分析
write.table(diff,file="edgerOut.xls",sep="\t",quote=F)
diffSig = diff[(diff$FDR < padj & (diff$logFC>foldChange | diff$logFC<(-foldChange))),]
write.table(diffSig, file="diffSig.xls",sep="\t",quote=F)
diffUp = diff[(diff$FDR < padj & (diff$logFC>foldChange)),]
write.table(diffUp, file="up.xls",sep="\t",quote=F)
diffDown = diff[(diff$FDR < padj & (diff$logFC<(-foldChange))),]
write.table(diffDown, file="down.xls",sep="\t",quote=F)
normalizeExp=rbind(id=colnames(newData),newData)
write.table(normalizeExp,file="normalizeExp.txt",sep="\t",quote=F,col.names=F) #输出所有基因校正后的表达值(normalizeExp.txt)
diffExp=rbind(id=colnames(newData),newData[rownames(diffSig),])
write.table(diffExp,file="diffmRNAExp.txt",sep="\t",quote=F,col.names=F) #输出差异基因校正后的表达值(diffmRNAExp.txt)
heatmapData <- newData[rownames(diffSig),]
#volcano火山图
pdf(file="vol.pdf")
xMax=max(-log10(allDiff$FDR))+1 #这一步很关键,如果allDiff里面有0,那xMax是无穷大,所以会报错,让你输入有限数值的定义,这里我们可以看一下我们的数据里面,FDR的最小值是多少,计算出来,直接输入该值,即xMax=300,就像下面的yMax。
yMax=12
plot(-log10(allDiff$FDR), allDiff$logFC, xlab="-log10(FDR)",ylab="logFC",
main="Volcano", xlim=c(0,xMax),ylim=c(-yMax,yMax),yaxs="i",pch=20, cex=0.4)
diffSub=allDiff[allDiff$FDR<padj & allDiff$logFC>foldChange,]
points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="red",cex=0.4)
diffSub=allDiff[allDiff$FDR<padj & allDiff$logFC<(-foldChange),]
points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="green",cex=0.4)
abline(h=0,lty=2,lwd=3)
dev.off()
#heatmap热图
hmExp=log10(heatmapData+0.001)
library('gplots')
hmMat=as.matrix(hmExp)
pdf(file="heatmap.pdf",width=60,height=90)
par(oma=c(10,3,3,7))
heatmap.2(hmMat,col='greenred',trace="none")
dev.off()
结果如下,我们产生了相关文件,这些文件都是在后续画图中有用的文件。

热图和火山图如下


下面我们进一步下载miRNA的数据,与上述不同的是,这次我们选择miRNAseq,而不是mRNAseq。分析完了以后看看关联情况。与上述下载方式类似。
把所有需要的文件下载下来,进一步进行分析。需要用到menifest文件。

cart文件需要代码下载。如上图所示
那么下载下来以后,我们继续把miRNA的文件融合到一起,这里不需要放到同一个文件夹,因为现在已经是txt文件,直接用perl脚本整合到一个文件夹。整合完了以后,已经是标准的miRNA ID,不再需要转换。下面进行差异分析和热图火山图的绘制。融合的文件如下,miRNAmatrix.txt,内容如下图。



上图显示,normal52,后面分组需要用到,请记录。
foldChange=1
padj=0.05
library("edgeR")
rt=read.table("miRNAmatrix.txt",sep="\t",header=T,check.names=F) #改成自己的文件名
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data=data[rowMeans(data)>0,]
#group=c("normal","tumor","tumor","normal","tumor")
group=c(rep("normal",4),rep("tumor",179)) #按照癌症和正常样品数目修改
design <- model.matrix(~group)
y <- DGEList(counts=data,group=group)
y <- calcNormFactors(y)
y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)
et <- exactTest(y,pair = c("normal","tumor"))
topTags(et)
ordered_tags <- topTags(et, n=100000)
allDiff=ordered_tags$table
allDiff=allDiff[is.na(allDiff$FDR)==FALSE,]
diff=allDiff
newData=y$pseudo.counts
write.table(diff,file="edgerOut.xls",sep="\t",quote=F)
diffSig = diff[(diff$FDR < padj & (diff$logFC>foldChange | diff$logFC<(-foldChange))),]
write.table(diffSig, file="diffSig.xls",sep="\t",quote=F)
diffUp = diff[(diff$FDR < padj & (diff$logFC>foldChange)),]
write.table(diffUp, file="up.xls",sep="\t",quote=F)
diffDown = diff[(diff$FDR < padj & (diff$logFC<(-foldChange))),]
write.table(diffDown, file="down.xls",sep="\t",quote=F)
normalizeExp=rbind(id=colnames(newData),newData)
write.table(normalizeExp,file="normalizeExp.txt",sep="\t",quote=F,col.names=F) #输出所有miRNA校正后的表达值(normalizeExp.txt)
diffExp=rbind(id=colnames(newData),newData[rownames(diffSig),])
write.table(diffExp,file="diffmiRNAExp.txt",sep="\t",quote=F,col.names=F) #输出差异miRNA校正后的表达值(diffmiRNAExp.txt)
heatmapData <- newData[rownames(diffSig),]
#volcano
pdf(file="vol.pdf")
xMax=max(-log10(allDiff$FDR))+1
yMax=12
plot(-log10(allDiff$FDR), allDiff$logFC, xlab="-log10(FDR)",ylab="logFC",
main="Volcano", xlim=c(0,xMax),ylim=c(-yMax,yMax),yaxs="i",pch=20, cex=0.4)
diffSub=allDiff[allDiff$FDR<padj & allDiff$logFC>foldChange,]
points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="red",cex=0.4)
diffSub=allDiff[allDiff$FDR<padj & allDiff$logFC<(-foldChange),]
points(-log10(diffSub$FDR), diffSub$logFC, pch=20, col="green",cex=0.4)
abline(h=0,lty=2,lwd=3)
dev.off()
#heatmap
hmExp=log10(heatmapData+0.001)
library('gplots')
hmMat=as.matrix(hmExp)
pdf(file="heatmap.pdf",width=60,height=30)
par(oma=c(10,3,3,7))
heatmap.2(hmMat,col='greenred',trace="none")
dev.off()


上图是miRNA的差异结果图。
下面,我们要进一步分析差异的mRNA和差异的miRNA有什么联系呢?
我们把diffmiRNAExp.txt和 diffmRNAExp.txt 放到目录文件夹,用perl脚本将其融合到一个文件。
miRNA = read.table("diffmiRNAExp.txt", row.names=1 ,header=T,sep="\t",check.names=F)
RNA = read.table("diffmRNAExp.txt", row.names=1 ,header=T,sep="\t",check.names=F)
colnames(miRNA)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*","\\1\\-\\2\\-\\3\\-\\4",colnames(miRNA))
colnames(RNA)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*","\\1\\-\\2\\-\\3\\-\\4",colnames(RNA))
sameSample=intersect(colnames(miRNA),colnames(RNA))
merge=rbind(id=sameSample,miRNA[,sameSample],RNA[,sameSample])
write.table(merge,file="merge.txt",sep="\t",quote=F,col.names=F)
miRNALabel=cbind(rownames(miRNA),"miRNA")
RNALabel=cbind(rownames(RNA),"gene")
nodeLabel=rbind(c("ID","Classify"),miRNALabel,RNALabel)
write.table(nodeLabel,file="nodeLabel.txt",sep="\t",quote=F,col.names=F,row.names=F)

nodelabel就是指定哪些属于miRNA,哪些属于mRNA,给一个标签,后续用cytoscape画图时候用,保存好。
然后下一步我们用WGCNA去分析共表达情况。
待续。。。。。。
请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!