转录组差异分析
差异分析难点在于将数据处理成需要的格式
表达矩阵
数值型矩阵-count
行名是symbol
低表达量的基因需要过滤
分组信息
因子，对照组在level第一位
与表达矩阵的列一一对应
项目名称
字符串（不要有特殊字符）
TCGA-XXX
非TCGA数据特殊无要求
-
拿不到count数据如何做差异分析？
-
自行做上游分析得到count
-
tpm：取log，用limma做差异分析
-
fpkm、rpkm：转换为tpm，取log，用limma做差异分析
-
- 不同类型转录组数据的应用
差异分析数据整理
差异分析的前提：count数据
项目取名字
TCGA的数据，统一叫TCGA-xxxx，非TCGA的数据随意起名，不要有特殊字符即可。
proj = "TCGA-CHOL"
表达矩阵
dat = read.table("TCGA-CHOL.htseq_counts.tsv.gz",check.names = F,row.names = 1,header = T)
range(dat)
#> [1] 0.0000 24.1811
#逆转log,发现需要逆转，才逆转
dat = as.matrix(2^dat - 1)
dat[1:4,1:4]
#> TCGA-ZD-A8I3-01A TCGA-W5-AA2U-11A TCGA-W5-AA30-01A
#> ENSG00000000003.13 5254 2476 5132
#> ENSG00000000005.5 1 1 0
#> ENSG00000000419.11 1212 655 1644
#> ENSG00000000457.12 753 346 2652
#> TCGA-W5-AA38-01A
#> ENSG00000000003.13 8249
#> ENSG00000000005.5 1
#> ENSG00000000419.11 1696
#> ENSG00000000457.12 519
# 深坑一个
dat[97,9]
#> [1] 876
as.character(dat[97,9]) #眼见不一定为实吧。
#> [1] "875.999999999999"# 转换为整数矩阵
exp = round(dat)
# 检查
as.character(exp[97,9])
#> [1] "876"
临床信息
clinical = read.delim("TCGA-CHOL.GDC_phenotype.tsv.gz")
clinical[1:4,1:4]
#> submitter_id.samples age_at_initial_pathologic_diagnosis
#> 1 TCGA-ZH-A8Y2-01A 59
#> 2 TCGA-ZH-A8Y7-01A 59
#> 3 TCGA-W7-A93O-01A NA
#> 4 TCGA-W7-A93O-11A NA
#> albumin_result_lower_limit albumin_result_specified_value
#> 1 NA NA
#> 2 3.5 2.4
#> 3 NA NA
#> 4 NA NA
表达矩阵行名ID转换
library(tinyarray)
exp = trans_exp_new(exp)
#> Warning in AnnoProbe::annoGene(rownames(exp), ID_type = "ENSEMBL", species =
#> species): 6.54% of input IDs are fail to annotate...
exp[1:4,1:4]
#> TCGA-ZD-A8I3-01A TCGA-W5-AA2U-11A TCGA-W5-AA30-01A TCGA-W5-AA38-01A
#> DDX11L1 0 0 0 1
#> WASH7P 81 10 146 55
#> MIR6859-1 1 0 11 1
#> MIR1302-2HG 0 0 0 0
基因过滤
需要过滤一下那些在很多样本里表达量都为0或者表达量很低的基因。过滤标准不唯一。
过滤之前基因数量：
nrow(exp)
#> [1] 56514
- 常用过滤标准1
仅去除在所有样本里表达量都为零的基因
exp1 = exp[rowSums(exp)>0,]
nrow(exp1)
#> [1] 48057
- 常用过滤标准2
仅保留在一半以上样本里表达的基因
exp = exp[apply(exp, 1, function(x) sum(x > 0) > 0.5*ncol(exp)), ]
nrow(exp)
#> [1] 28434
分组信息获取
TCGA的数据，直接用make_tcga_group给样本分组（tumor和normal），其他地方的数据分组方式参考芯片数据pipeline/02_group_ids.R
library(tinyarray)
Group = make_tcga_group(exp)
table(Group)
#> Group
#> normal tumor
#> 9 36
保存数据
save(exp,Group,proj,clinical,file = paste0(proj,".Rdata"))
玩转GEO的实用工具
library(tinyarray)
get_count_txt('GSE204753')
##获得超级标准的表达矩阵！！
引用自生信技能树课程，又是爱小洁老师的一天！！
