基因课FTP地址：ftp://http://gsx.genek.tv/2020-3-10%E7%9B%B4%E6%92%AD%E4%B8%80%E4%B8%AA%E5%AE%8C%E6%95%B4%E7%9A%84%E8%BD%AC%E5%BD%95%E7%BB%84%E9%A1%B9%E7%9B%AE/
听张旭东老师的课

加载tidyverse包

library(tidyverse)

环境数据保存与加载

保存
save(gene_exp, sample_info, gene_info, file = 'data/rnaseq-apple/rprepare.rdata')
加载
load(file = 'data/rnaseq-apple/rprepare.rdata')

导入差异表达分析结果

de <- read.table(file = 'data/rnaseq-apple/genes.counts.matrix.KID_S1_vs_KID_S4.DESeq2.DE_results', header = T) # 不将第一列设置为列名，因为要用tidyverse处理，tidyverse体系中不需要行名

提取想要的列(注意：有些绘图要求不去除小于阈值的FC值)

查看所有列名
colnames(de)
选择列
select(de, id, log2FoldChange, pvalue, padj) 或
select(de, -baseMeanA, -baseMeanB, -baseMean, -lfcSE, -stat) # 列名前加 ‘-’减号表示去除该列
过滤行
filter(de, abs(log2FoldChange) > 1 & padj < 0.05)
两个操作同步完成
test <- select(de, id, log2FoldChange, pvalue, padj) %>%
filter(abs(log2FoldChange) > 1 & padj < 0.05)

新建列表示FC(表达量倍数)、direction(表达量升高还是降低)

deg <- mutate(FC = 2 ** log2FoldChange, direction = if_else(log2FoldChange > 1, 'up', 'down'))

多表关联聚合LEFT JOIN

Excel中用lookup
R语言中
left join(deg, gene_info, by = c('id' = 'GID')) %>% # 关联基因信息
left join(rownames_to_column(gene_exp, var = 'id'), by = 'id') # 关联表达量信息

整合后的步骤

de_result <- # 数据导入
mutate(de_result, direction =
if_else(padj > 0.05, 'ns',
if_else(abs(log2FoldChange) < 1, 'ns', # 同时满足padj<0.05, log2FC>1才能叫significant
if_else(log2FoldChange >= 1, 'up', 'down'))) # 添加上下调信息
) %>%
left join(gene_info, by = c('id' = 'GID')) %>% # 关联基因信息
left join(rownames_to_column(gene_exp, var = 'id'), by = 'id') %>% # 关联表达量信息
dplyr::select(-c(2:4, 6:7)) %>% # 去除无用的列
arrange(desc(abs(log2FoldChange))) # 按log2FoldChange绝对值降序排列

按照上调与下调分组

统计上调、下调、ns的基因数量
group_by(de_result, direction) %>%
summarise(count = n())
找到上调下调最大的基因
group_by(df, ID) %>%
filter(FC == max(FC))

题外话

用vim 做文件中内容替换
: → s/要替换的字符/替换后的字符/ # 替换光标所在行的第一个匹配的要替换的字符
: → s/要替换的字符/替换后的字符/g # 替换光标所在行中所有的对应字符
: → %/要替换的字符/替换后的字符/g # 替换文本中所有的对应字符

Tidyverse处理差异分析结果