起源是海洋师兄让我复现一个公司做的生信分析，据说这样一个东西收费1w+，于是我开始尝试进行探索。
分析GEO数据miRNA表达情况
GSE106452是4种肝癌细胞的外泌体miRNA的微阵列分析数据，我们由于只有HepG2细胞，因此只拿该细胞数据，挑选表达较高的30个miRNA数据。甚至都用不到Python和R，Excel排序即可。
点评：开始俺还以为得用GEO2R做差异分析，事实上是俺想复杂了。
TCGA差异分析
TCGA数据下载
选择数据
访问GDC，在Cases里的Project选择TCGA-LIHC，Files里的Data Category选择transcriptome profiling，如果需要下载临床资料，类似地，Data Category选择Clinical即可。
原始方法
如果大小不超过50M，可以点击右边的Add All Files to Cart，然后再右上角的Cart里面选择Download即可下载。
下载GDC Data Transfer Tool，并配置环境变量。
在之前的页面Data Category选择transcriptome profiling，Experimental Strategy选择RNA-Seq，然后在右边选择下载Manifest，文件命名为gdc_manifest.2020-01-17-LIHC-RNA-Seq.txt，此外Experimental Strategy选择miRNA-Seq，Data Type分别选择miRNA Expression Quantification和Isoform Expression Quantification分别得到另外两个文件gdc_manifest.2020-01-17-LIHC-miRNA-Seq.txt和gdc_manifest.2020-01-17-LIHC-miRNA-Isoform.txt。类似地，Data Category选择Clinical，Data Format选择bcr xml即可获取gdc_manifest.2020-01-17-LIHC-Clinical.txt。
新建一个LIHC文件夹，将上述Manifest拷进去，CMD下输入gdc-client download -m gdc_manifest.2020-01-17-LIHC-Clinical.txt -d Clinical/ --log-file Clinical.log即可下载临床数据，其他的操作同理。不过无力吐槽这个网络可连接性，下载起来比登天还难。
这里我采用了分割文件&挂代理的策略才把RNA数据下完。
console下设置代理：
set http_proxy=http://127.0.0.1:10809。
set https_proxy=https://127.0.0.1:10809。
此外，在下载上述manifest文件以外顺便应该下载json文件备用。
数据处理
Clinical数据合并
主要是xml文件，利用R包来分析。以单个文件为例：
1
2
3
4
5
6
7
| library("XML")
library("methods")
result = xmlParse(file="D:/Bioinformatics/TCGA/LIHC/Clinical/00a9a7f4-06eb-40fb-8d3a-66f5f5d315f7/nationwidechildrens.org_clinical.TCGA-KR-A7K0.xml")
rootnode = xmlRoot(result)
rootsize = xmlSize(rootnode)
xmldataframe = xmlToDataFrame(rootnode[2])
write.table(t(xmldataframe),'tmp')
|
当然具体数据有成百上千个，应该用循环来做。
1
2
3
4
5
6
7
8
9
10
11
| library("XML")
library("methods")
dir = "D:/Bioinformatics/TCGA/LIHC/Clinical/"
cl = lapply(list.files(path=dir,pattern="*.xml$",recursive=T),function(x){
result = xmlParse(file=file.path(dir,x))
rootnode = xmlRoot(result)
xmldataframe = xmlToDataFrame(rootnode[2])
return(t(xmldataframe))
})
cl_df = t(do.call(cbind,cl))
save(cl_df,file="D:/Bioinformatics/TCGA/LIHC/Process/GDC_TCGA_LIHC_clinical_df.Rdata")
|
miRNA-Seq数据合并
类似地，修改处理代码如下：
1
2
3
4
5
6
7
8
9
| dir = "D:/Bioinformatics/TCGA/LIHC/miRNA-Seq/"
mi = lapply(list.files(path=dir,pattern="*.mirnas.quantification.txt$",recursive=T),function(x){
result = read.table(file=file.path(dir,x),sep="\t",header=T)[,1:2]
return(result)
})
mi_df = t(do.call(cbind,mi))
colnames(mi_df) = mi_df[1,]
mi_df = mi_df[seq(2,nrow(mi_df),by=2),]
save(mi_df,file="D:/Bioinformatics/TCGA/LIHC/Process/GDC_TCGA_LIHC_miRNA-Seq_df.Rdata")
|
这样可以得到一个表达矩阵。
Json数据解析
类似地，修改处理代码如下：
1
2
3
4
5
6
7
8
9
| library("rjson")
result = fromJSON(file = "D:/Bioinformatics/TCGA/LIHC/TCGA-LIHC-miRNA-Seq.json")
fls = unlist(lapply(result,function(x){x$file_name}))
cid = unlist(lapply(result,function(x){x$cases[[1]]$case_id}))
id2fls = data.frame(cid=cid,fls=fls)
save(id2fls,mi_df,fls,cl_df,file='D:/Bioinformatics/TCGA/LIHC/Rdata/GDC_TCGA_LIHC_miRNA-clinical.Rdata')
rm(list=ls())
load(file='D:/Bioinformatics/TCGA/LIHC/Rdata/GDC_TCGA_LIHC_miRNA-clinical.Rdata')
|
使用TCGAbiolinks包
感觉真的很方便，在BiocManager里面安装：
1
2
3
4
5
6
7
8
9
10
11
12
|
BiocManager::install("TCGAbiolinks")
library(TCGAbiolinks)
query = GDCquery(project = 'TCGA-LIHC',
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - Counts")
GDCdownload(query, method = "api", files.per.chunk = 100)
expdat = GDCprepare(query = query)
count_matrix = assay(expdat)
|
PS1：本来打算用RTCGA的，但是这个包N年不更新了，用着膈应~
PS2：R语言设置代理：
Sys.setenv(“http_proxy”=”http://127.0.0.1:10809")
Sys.setenv(“https_proxy”=”https://127.0.0.1:10809")
##
