01072019 更新
所以说单位发了一个顶配的Mac Pro，然后随意Run了下，轻松读进去了。
rt<-read.table("GTEx_Analysis_2016-01-15_v7_RNASeQCv1.1.8_gene_tpm.gct", skip = 2, header = TRUE, sep = "\t")
=============
GTEx数据库不用多说，我下载下来了表达数据为gct格式文件但是特别大，怎么办。脚本语言perl和python都刚学还不会，所以就准备用R来操作。
我的目的非常简单，提取GTEx中肝组织的表达数据。因此我下载了gct的表达数据，和组织的生物学信息从中获取了liver组织的ID。下载网址为：https://gtexportal.org/home/datasets
要注意代码中，ID编号的连接符的一致性因为R会默认把列名中的连接符-变为.,所以要注意替换 或者要加上check.names=F。
我使用了两个文件：liver.csv这里是肝细胞的ID，从GTEx下载的数据提取出来的，GTEx_Analysis_2016-01-15_v7_RNASeQCv1.1.8_gene_tpm.gct这个是GTEx下载的表达数据。
我参考了这个教程：https://blog.csdn.net/u012432611/article/details/50224015
然后自己写了个循环构建了classes变量。
liverid<-read.csv("liver.csv",header=T)
liverid<-as.character(liverid[,1])
rt<-read.table("GTEx_Analysis_2016-01-15_v7_RNASeQCv1.1.8_gene_tpm.gct", skip = 2, header = TRUE, sep = "\t",nrow=5,check.names=F)
liverid<-intersect(colnames(rt),liverid)
c<-vector()
a<-colnames(rt)
for (i in 1:length(liverid)){
b<-which(a==liverid[i])
c<-c(c,b)
}
d<-c
classes <- sapply(rt, class)
classes[d] <- rep("NULL", length(classes)-length(d))
rtt<-read.table("GTEx_Analysis_2016-01-15_v7_RNASeQCv1.1.8_gene_tpm.gct", skip = 2, header = TRUE, sep = "\t",colClasses=classes, check.names=F)
write.csv(rtt,"liver_GTEx_expression.csv")
网友评论