基因课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 # 替换文本中所有的对应字符
网友评论