生信马拉松 Day21 转录组的分析实战

啊啊啊，太伤心了，这一天的课小洁老师抽了我的数据集做师范，我竟然上一半跑路么有上和甜甜的小洁连麦的机会o(╥﹏╥)o

今天主要是实战演练，顺便复习了R的函数以及Rmarkdown的用法

内容一：R函数的复习

这里分享我不太熟练的stringr包的函数

# 编一个字符串
library(stringr)
a = c("A_con_1","A_con_2","A_con_3","A_ZY_1","A_ZY_2","A_ZY_3","ES2con_1","ES2con_2","ES2con_3","ES2ZY_1","ES2ZY_2","ES2ZY_3")
#删除最后两个字符
library(stringr)
#正则表达式，\\d
str_remove(a,'_\\d')
str_sub(a,1,-3)
#字符串替换
str_replace(a,"ES2","ES2_")
#检测是否有ZY字符
str_detect(a,'ZY')

内容二：转录组数据

首先需要下载自己对应的数据，GEO的转录组数据标识为Expression profiling by high throughput sequencing，需要注意单细胞转录组也被纳入其中，需要自己看详情页区分，单细胞转录组不能用下面的数据处理方法

其次，我们做转录组差异分析用的是count值，这可以在样本详情页寻找对数据的注释信息，或者下载Supplementary file文件解压打开之后是整数（除非有对数据的特别解释说明）

注意不能照搬前面芯片分析的过程，因为转录组和芯片的差异的技术手段和来源不一样，数据的含义有差别，所以处理也不同

count/reads计数数据

只有转录组有count，芯片是表达量数据值

转录组数据在下机的时候，也就是从实验品变数据的时候需要软件来数，下机数据是短的序列片段reads形如ATCG等，通过比对给各个短片段找到位置，一个基因的管辖范围内找到多少reads就是多少，每个基因和样本都得到这个count数就形成矩阵，因此可能有基因没有对应的片段就是0

来自：生信技能树，生信马拉松，小洁老师

count是转录组结果的格式之一，没有count数值没法做差异分析，tpm、fpkm、rpkm都是对count值的转换，这些无法转化回count，最好的选择是原始的count

没有count值的处理方法：

1.自行做上游分析得到count矩阵

2.tpm：取log，用limma做差异分析（迫不得已）

3.fpkm、rpkm：转换为tpm，取log，用limma做差异分析（迫不得已）

可参考：https://mp.weixin.qq.com/s/_DtkxSfLGQHcRju66J4yTQ

另外关于expected_count和norm_count，注意即edgeR只能用expected，vomm理论上可以使用norm_count（只是可以不是必须）。参考https://www.jianshu.com/p/46b048220b88

来自：生信技能树，生信马拉松，小洁老师

转录组输入数据的要求，来自：生信技能树，生信马拉松，小洁老师

转录组的输入数据是来自补充文件里，内容格式不确定，目标是变成count矩阵，行名是基因名称，列名只要是不同的就行。整理的过程比较困难，不像芯片有exprs可以直接提取

差异分析有3个包进行差异分析

DESeq2

edgeR

limma

三个包都值得学习，虽然名字和函数不同，结果都是logFC和p.value

三个包都在count的基础上进行标准化处理，然后进行logFC转化，所以3个包的差异基因不同

三个R包就会有3组差异基因，用韦恩图展示交集

cpm，tpm，fpkm，rpkm都是log后用的，可以进行pca、生存分析、模型构建、聚类分析、相关性等，但除了差异分析

整理数据的主要矛盾：只需认准count数据即可，自己的数据还是公共数据、数据库、背景知识均不影响差异分析，xxxm的数据需要转化为tpm后走芯片分析的流程

内容三：示例数据分析流程

1.起个项目名字

TCGA的数据，统一叫TCGA-xxxx，非TCGA的数据随意起名，不要有特殊字符即可。

proj = "TCGA-CHOL"  

2.读取和整理数据

2.1 表达矩阵

dat = read.table("TCGA-CHOL.htseq_counts.tsv.gz",check.names = F,row.names = 1,header = T)
#注意这里的参数酌情选择，特别是作为行名的列有重复以及列名中有特殊字符时
range(dat)
#取过log的数据一般在20以内，正常的数据几十几百几千都有且是整数
#需要注意R在显示的时候会在一定位数之后近似，使0.9999999显示出来是1
#这时可以用as.charcter()查看，此时数字可以看到数据真实的样子
#本例的存在两个问题：
#1.数据看网页介绍做过log，需要逆转
#2.dat的行名不是symbol需要转化
dat = as.matrix(2^dat - 1)
dat[1:4,1:4]
# 深坑一个
dat[97,9]
as.character(dat[97,9]) #眼见不一定为实吧。

# 转换为整数矩阵
exp = round(dat)
# 检查
as.character(exp[97,9])

2.2 临床信息

clinical = read.delim("TCGA-CHOL.GDC_phenotype.tsv.gz")
clinical[1:4,1:4]

3.表达矩阵行名ID转换

library(tinyarray)
exp = trans_exp_new(exp) 
#如果这个函数不得行，按住ctrl再点这个函数可以看到源代码，然后自己根据自己的数据操作
#有一些fail正常
#默认物种是人类，需要修改物种选参数species
exp[1:4,1:4]

4.基因过滤

需要过滤一下那些在很多样本里表达量都为0或者表达量很低的基因。过滤标准不唯一。

过滤之前基因数量：

nrow(exp)

常用过滤标准1：

仅去除在所有样本里表达量都为零的基因

exp1 = exp[rowSums(exp)>0,]
nrow(exp1)

常用过滤标准2(推荐)：

仅保留在一半以上样本里表达的基因

exp = exp[apply(exp, 1, function(x) sum(x > 0) > 0.5*ncol(exp)), ]
#这一行理解有个难点，这儿不是sum(x)即对行求和而是sum(x>0)，类似的sum(32>0)的返回值是1
nrow(exp)

5.分组信息获取

TCGA的数据，直接用make_tcga_group给样本分组（tumor和normal），其他地方的数据分组方式参考芯片数据pipeline/02_group_ids.R

library(tinyarray)
Group = make_tcga_group(exp)
table(Group)

6.保存数据

save(exp,Group,proj,clinical,file = paste0(proj,".Rdata"))

第二部分的代码只需要修改输入数据的名称即可

内容四：当GEO上数据不是count也转不回count的时候

NCBI整理的数据没有正常的表达矩阵，但是又不想搞上游分析时候的方法

library(tinyarray)
get_count_txt("GSE190518")
#会回复一个网页，把网页复制到浏览器里就可以看
#目前只有人类的可以

注意这个写进R markdown文件里时，若设置了knitr的message=F就看不到了

这种方法出来的样本数可能和原始的丢失样本，因为NCBI会对样本质量进行过滤，是正常的

生信技能树，生信马拉松，小洁老师