已测试：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)