复现文献：doi: 10.7150/ijbs.41587
gdc数据整理
1.表达矩阵整理
step1.提取数据
#step0.加载R包
library(tibble)
library(dplyr)
library(stringr)
#step1.读取metadata数据----
rm(list = ls())
metaDat <- jsonlite::fromJSON("../01_exprCount/metadata.cart.2020-09-21.json")
#step2.获取样本名与文件名的对应关系----
##file_id是文件夹的名字
##file_name是文件夹压缩包的名字,即表达数据压缩包
##associated_entities里面包含对应的样本名
filesPath <- paste("../01_exprCount",metaDat$file_id,
metaDat$file_name,sep = "/") #创建表达文件路径,../当前目录上一级
filesPath[1] #查看是否正确
# "../01_exprCount/80f895da-7635-47a4-b974-93218702c9a5/7f3b65a3-ad59-49f1-92fa-72793435f86e.htseq.counts.gz"
tmp <- data.table::fread(filesPath[1]);head(tmp) #测试是否成功读取数据
rm(tmp) #节省空间，耳目一新
class(metaDat$associated_entities) #数据结构是list
sampleName <- metaDat$associated_entities%>%
unlist()%>% #破坏list结构
.[grep("^TCGA.*",.)]%>% #获取样本名
as.character();length(sampleName) #将样本名变为字符串,确定字符串长度
metaDat$associated_entities[[1]];sampleName[1] #查看顺序是否改变
#step3.测试，准备批量处理----
test1 <- data.table::fread(filesPath[1]);head(test1)
test2 <- data.table::fread(filesPath[2]);head(test2)
#identical(test1$V1,test2$V1) #可以鉴定，也可以不鉴定，因为后面是merge
test <- merge(test1,test2,by="V1");head(test)
rm(test,test1,test2) #节省空间，耳目一新
#step4.批量读取整理数据----
n=0
for(i in filesPath){
n=n+1
text = data.table::fread(i)
colnames(text)[2]=sampleName[n]
if(n<2){
exprDat=text
}else{
exprDat=merge(exprDat,text,by="V1")
}
}
exprDat[1:4,1:4] #查看是否成功提取
colnames(exprDat)[1] <- "ensemblID" #整理列名
str_remove_all("ENSG00000000003.13","\\..*") #不确定时，先测试一下
exprDat$ensemblID <- str_remove_all(exprDat$ensemblID,"\\..*");exprDat$ensemblID[1:3]
exprDat <- column_to_rownames(exprDat,var = "ensemblID")
save(exprDat,file = "exprDat_read.Rdata") #对于处理时间较长的数据都建议先保存，避免下面出bug
一定要验证！！！
小坑：误以为以manifest为基准进行Merge，得到的数据与manifest数据排序一样，不会改变
filenames <- paste(manifest$id,manifest$filename,sep = "/")
错误结果：
step2.滤过低表达基因
#step5.去除低表达量基因----
count0 <- apply(exprDat, 1, function(x)sum(x==0)) #统计count为0的情况
hist(count0,breaks = 100,col = "#efc000",main = "count为0分布情况") #可视化
dim(exprDat) #60488 373
lowLevel <- floor(ncol(exprDat)*0.1) #确定低表达的阈值，ceiling也可以
exprDat <- exprDat[apply(exprDat, 1, function(x) sum(x==0) < lowLevel ), ]
dim(exprDat) #26019 373
Q1:为什么要滤过低表达基因
A1：
-
从生物学的角度来看，在任何情况下都没有以生物学意义表达的基因不会引起人们的关注，因此最好将其忽略掉，因为这会增加总数 经过多次测试校正后的DEG，提高了过滤后DEG的灵敏度和精度。
-
从数据角度来看，通常将读取计数在0-10范围内的基因/转录物视为人工产物或“噪音”。
-
从统计学的角度来看，删除低计数基因可以更可靠地估计数据中的均值-方差关系。
-
从操作角度来看，删除低表达量的基因，有利于差异分析时，减少构建dss对象的时间。
Q2:如何滤过低表达基因
A2：见参考
参考：https://seqqc.wordpress.com/2020/02/17/removing-low-count-genes-for-rna-seq-downstream-analysis/
step3.创建分组信息
#step6.创建分组信息----
groupList <- ifelse(as.numeric(str_sub(colnames(exprDat),14,15))<10,
"tumor","normal")
groupList <- factor(groupList,levels = c("normal","tumor"))
groupList
table(groupList)
[options]分离lncRNA&mRNA
为什么需要分离lncRNA和mRNA呢？
这个是根据不同的分析需求决定的需不需要分离，比如构建ceRNA网络时，需要进行分离操作！
step1.数据准备:count表达矩阵&注释文件(gtf/TCGA_anno.Rdata)
下载注释文件，最好是官网指定的注释文件（也可以在GENCODE下载）;注意lncRNA的定义
step2.分离数据:获取注释文件&分离
#step7.分离lncRNA & mRNA----
##step7.1读取并探索gtf文件
options(stringsAsFactors = F)
library(rtracklayer)
gtfDat <- rtracklayer::import("gencode.v22.annotation.gtf") #燃烧小电脑
class(gtfDat)
# [1] "GRanges"
# attr(,"package")
# [1] "GenomicRanges"
gtfDat <- as.data.frame(gtfDat);dim(gtfDat)
# [1] 2563671 27
colnames(gtfDat)
# [1] "seqnames" "start"
# [3] "end" "width"
# [5] "strand" "source"
# [7] "type" "score"
# [9] "phase" "gene_id"
# [11] "gene_type" "gene_status"
# [13] "gene_name" "level"
# [15] "havana_gene" "transcript_id"
# [17] "transcript_type" "transcript_status"
# [19] "transcript_name" "tag"
# [21] "transcript_support_level" "havana_transcript"
# [23] "exon_number" "exon_id"
# [25] "ont" "protein_id"
# [27] "ccdsid"
fetureType <- as.data.frame(table(gtfDat$type));fetureType
# Var1 Freq
# 1 gene 60483
# 2 transcript 198442
# 3 exon 1172082
# 4 CDS 699443
# 5 start_codon 82228
# 6 stop_codon 74337
# 7 UTR 276542
# 8 Selenocysteine 114
geneType_all <- as.data.frame(table(gtfDat$gene_type));geneType_all
##step7.2:先筛选出gene对应的行
nrow(gtfDat);gtfDat <- gtfDat[gtfDat$type=="gene",];nrow(gtfDat)
geneType <- as.data.frame(table(gtfDat$gene_type));geneType#查看排除非基因组件后的基因类型
##step7.3:提取lnc和mRNA及其对应的ensambelid和symbol
lncType <- c("3prime_overlapping_ncRNA", "antisense",
"bidirectional_promoter_lncRNA", "lincRNA",
"macro_lncRNA", "non_coding",
"processed_transcript", "sense_intronic" ,
"sense_overlapping")
#一定要注意lncRNA的范围：https://www.gencodegenes.org/pages/biotypes.html
table(gtfDat$gene_type %in% lncType) #探索lncRNA数量
table(gtfDat$gene_type == "protein_coding") ##探索mRNA数量
lncAnno <- gtfDat[gtfDat$gene_type %in% lncType,c("gene_name","gene_id","gene_type")]
mAnno <- gtfDat[gtfDat$gene_type == "protein_coding",c("gene_name","gene_id","gene_type")]
head(lncAnno);head(mAnno)
lncAnno$gene_id <- str_remove_all(lncAnno$gene_id,"\\..*")
mAnno$gene_id <- str_remove_all(mAnno$gene_id,"\\..*")
head(lncAnno);head(mAnno)
save(lncAnno,mAnno,file = "lncRNA_mRNA_anno.Rdata")
#保存好注释文件后，在官方未作出修改说明的前提下，以后的分析不需要再读取gtf文件，直接加载Rdata数据即可
##step7.4:表达矩阵拆分和注释
#初步探索一下
exprDat <- rownames_to_column(exprDat,var = "gene_id")
mRNA_exprDat <- merge(exprDat,mAnno,by="gene_id")
table(!duplicated(mRNA_exprDat$gene_name))
library(dplyr)
#开始
load("exprDat_read.Rdata")
mRNA_exprDat <- exprDat %>%
rownames_to_column(var = "gene_id") %>%
merge(mAnno,by="gene_id") %>% #合并探针的信息
select(-c("gene_id","gene_type")) %>% #去掉多余信息
select(gene_name,everything()) %>% #重新排列，
mutate(rowMean =rowMeans(.[grep("TCGA", names(.))])) %>% #求平均数
arrange(desc(rowMean)) %>% #把表达量的平均值按从大到小排序
distinct(gene_name,.keep_all = T) %>% # gene_name留下第一个
select(-rowMean) %>% #去除rowMean这一列
column_to_rownames(var = "gene_name")
mRNA_exprDat[1:4,1:4]
#lncRNA同mRNA
lnc_exprDat <- exprDat %>%
rownames_to_column(var = "gene_id") %>%
merge(lncAnno,by="gene_id") %>% #合并探针的信息
select(-c("gene_id","gene_type")) %>% #去掉多余信息
select(gene_name,everything()) %>% #重新排列，
mutate(rowMean =rowMeans(.[grep("TCGA", names(.))])) %>% #求出平均数
arrange(desc(rowMean)) %>% #把表达量的平均值按从大到小排序
distinct(gene_name,.keep_all = T) %>% # gene_name留下第一个
select(-rowMean) %>% #去除rowMean这一列
column_to_rownames(var = "gene_name")
lnc_exprDat[1:4,1:4]
save(mRNA_exprDat,lnc_exprDat,groupList,file = "lncRNA_mRNA_exprDat.Rdata")
2.整理临床数据
#step1.设置工作环境
getwd()
setwd("C:\\Users\\luobo\\Documents\\生信学习-R语言\\XPO1 pan_cancer\\ACC_TCGA\\clinical")
library(XML)
#step2.测试
result <- xmlParse("01e4b744-ba8e-488c-b026-ebf94a2ef977/nationwidechildrens.org_clinical.TCGA-OR-A5JJ.xml")
class(result)
# [1] "XMLInternalDocument" "XMLAbstractDocument"
rootnode <- xmlRoot(result)
rootsize <- xmlSize(rootnode);rootsize
# [1] 2 表示有两个节点
print(rootnode[1])#节点1，我们不需要
print(rootnode[2])
xmldataframe <- xmlToDataFrame(rootnode[2])
head(t(xmlToDataFrame(rootnode[2])))
# [,1]
# additional_studies ""
# tissue_source_site "OR"
# patient_id "A5JJ"
# bcr_patient_barcode "TCGA-OR-A5JJ"
# bcr_patient_uuid "9EC86C06-E7A9-4EA6-B36D-6F832909054B"
# informed_consent_verified "YES"
#step3.正式提取
xmls = dir(".",pattern = "*.xml$",recursive = T)#recursive表示是否递归目录
xmls[1]
# "01e4b744-ba8e-488c-b026-ebf94a2ef977/nationwidechildrens.org_clinical.TCGA-OR-A5JJ.xml"
transXmls = function(x){
result <- xmlParse(x)
rootnode <- xmlRoot(result)
xmldataframe <- xmlToDataFrame(rootnode[2])
return(t(xmldataframe))
}
cl = lapply(xmls,transXmls)
cl_df <- t(do.call(cbind,cl))
cl_df[1:3,1:3]
# additional_studies tissue_source_site patient_id
# [1,] "" "OR" "A5JJ"
# [2,] "" "OR" "A5LO"
# [3,] "" "OR" "A5K4"
clinical = data.frame(cl_df)
clinical[1:4,1:4]
# additional_studies tissue_source_site patient_id bcr_patient_barcode
# 1 OR A5JJ TCGA-OR-A5JJ
# 2 OR A5LO TCGA-OR-A5LO
# 3 OR A5K4 TCGA-OR-A5K4
# 4 OR A5JK TCGA-OR-A5JK
作者：骆栢维
编辑：骆栢维 周昕雨
校审：梁晓杰