所有的学习都是从模仿开始的,所以我又开始跟着小洁老师继续学习了。知识是日积月累的,希望自己走的每一步都是扎实的。
1.下载包及引用数据
#下载包和引用包
if (!require(RTCGA))BiocManager::install("RTCGA")
if (!require(RTCGA.miRNASeq))BiocManager::install("RTCGA.miRNASeq")
if (!require(RTCGA.clinical))BiocManager::install("RTCGA.clinical")
library(RTCGA.miRNASeq)
KIRC.miRNASeq[1:10,1:4]#使用KIPC的数据集
让我们看看他的结构长什么样,与别的方法对比一下是否一样

expr1 <- expressionsTCGA(KIRC.miRNASeq)#处理数据
dim(expr1)#计算维度
#[1] 1779 1049
expr1[1:10,1:4]#查看部分结构

结合 图KIRC.miRNASeq以及expr1[1:10,1:4]发现:
从第三列开始是miRNA表达量
每三行是一个样本的信息,而我们做数据只需要read_count
expressionsTCGA()函数处理后行名丢失,需要补充回去
2.整理为表达矩阵
根据上面发现的结果,取3至ncol()列,每三行取第一行,补上行名
expr2 = expr1[seq(1, nrow(expr1), by = 3), 3:ncol(expr1)]
dim(expr2)
#> [1] 593 1047
expr2[1:4,1:4]
#> # A tibble: 4 x 4
#> `hsa-let-7a-1` `hsa-let-7a-2` `hsa-let-7a-3` `hsa-let-7b`
#> <chr> <chr> <chr> <chr>
#> 1 5056 10323 5429 17908
#> 2 14503 29238 14738 37062
#> 3 8147 16325 8249 28984
#> 4 7138 14356 7002 6909
#此时表达量均为字符型,改为数值型
expr2 = apply(expr2, 2, as.numeric)
#Warning message:
#In apply(expr2, 2, as.numeric) : 强制改变过程中产生了NA
#行名补回去
rownames(expr2) = rownames(KIRC.miRNASeq)[seq(1, nrow(KIRC.miRNASeq), by = 3)]
expr = t(expr2)#转换行列
expr = na.omit(expr)#过滤NA
expr[1:4,1:4]

3.过滤低表达量的miRNA
过滤的标准:在10个以上样本中表达量>1,x>1返回逻辑值,sum()函数处理逻辑值向量,返回结果为TRUE的个数。
dim(expr)
#> [1] 1046 593
expr = expr[apply(expr, 1, function(x) {
sum(x > 1) > 10
}), ]
dim(expr)
#> [1] 552 593 #过滤掉454个miRNA
三.处理临床信息
library(RTCGA.clinical)
clinical <- KIRC.clinical
dim(clinical)
#> [1] 537 2809
#从2809列中跳出想要的列
clinical = clinical[c(
'patient.bcr_patient_barcode',
'patient.vital_status',
'patient.days_to_death',
'patient.days_to_last_followup',
'patient.race',
'patient.age_at_initial_pathologic_diagnosis',
'patient.gender' ,
'patient.stage_event.pathologic_stage'
)]
rownames(clinical) <- NULL
clinical <- tibble::column_to_rownames(clinical,var = 'patient.bcr_patient_barcode')#把指定的列转换为行名
rownames(clinical) <- stringr::str_to_upper(rownames(clinical))#实现大小写的转换
dim(clinical)
#> [1] 537 7
clinical[1:4,1:4]

![clinical[1:4,1:4]](https://img.haomeiwen.com/i21455872/
ee54529bb1d01537.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
保存数据
save(expr,file = "expr.Rda")
save(clinical,file = "clinical.Rda")
注:通过RTCGA包得到的数据并不是最新的,只是为了方便我学习和对比TCGA-biolinks包下载的数据。
结果:运用RTCGA包,我们获取到了miRNA表达矩阵和临床信息。表达矩阵中有593个样本,553个miRNA信息(过滤后的)。临床信息有537个样本,7个临床特征。
参考文章:TCGA-5.使用RTCGA包获取数据
网友评论