让我们来一起看看最热门的组学分析
各位解螺旋的小伙伴们大家周日好，欢迎大家来到每周日的先锋宇专栏！又是一周过去啦，放假的脚步越来越近，然后疫情防控的需要，相信很多小伙伴们也是只能待在当地过年吧。俗话说：年货可以囤，但是SCI不能屯。所以不能回家的日子就让先锋宇老师来陪伴大家继续学习科研知识，争取在过年之前赶紧弄出一篇SCI来投稿！
写在前面
通过前面几期的学习，我们已经掌握了基因组和转录组的分析方法，那从今天开始，先锋宇老师将带大家进入另外一个热门组学的分析：那便是免疫组学！毫无疑问，免疫肯定是当下最热门的组学，感觉无论做什么分析，都要和免疫挂上钩，这样才能给文章加点分。那接下来先锋宇老师就带大家看一看免疫组学到底有什么特别的地方让这么多人流连忘返，久久不忍离去。
代码演示
首先我们来复习一下前面介绍过的转录组数据的读取和清洗。
首先我们仍然要读取从Xena上面下载的胃癌数据到R里面：
library(tidyverse)stad_fpkm <-read_tsv("TCGA-STAD.htseq_fpkm.tsv.gz")

接下来我们指提取基因的表达量并且对表达矩阵进行基因的注释：
gtf_v22 <- read_tsv("gencode.gene.info.v22.tsv") %>% dplyr::filter(gene_type == "protein_coding")%>% dplyr::select(1,2)colnames(stad_fpkm)[1] <- "gene_id"stad_fpkm1 <- inner_join(gtf_v22,stad_fpkm, by = "gene_id") %>% dplyr::select(-1)

然后接下来我们把FPKM数据转换为TPM数据：
stad_fpkm1 <- aggregate(.~ gene_name, stad_fpkm1, mean)stad_fpkm1 <- stad_fpkm1 %>% column_to_rownames("gene_name")stad_fpkm1 <- 2^stad_fpkm1-1fpkmToTpm <- function(fpkm){ exp(log(fpkm) - log(sum(fpkm)) + log(1e6))}stad_tpm <- apply(stad_fpkm1, 2, fpkmToTpm)

然后接下来我们去除正常组织的样本，只保留肿瘤组织的样本：
data <- stad_tpm[,str_sub(colnames(stad_tpm),14,15) == "01"]

可以看到我们现在得到了一个19712行和375列的表达矩阵。
这些内容属于复习内容，所以讲得比较浅，想看更细致的讲解可以参考先锋宇专栏前面的推文。
免疫组学
那接下来我们就开始往免疫组学方向进行分析啦，首先我们对数据进行过滤，我们将基因在每个样本中都为0的那些基因给去掉，因为这样的基因没有任何意义。
data=data[rowMeans(data)>0,]

然后接下来我们使用limma包中的voom函数对我们的数据进行标准化，大概的意思也就是对我们的数据进行log2转化。
v <-voom(data,plot = F,save.plot = F)

然后我们从voom转化后的数据列表对象v中把我们的表达矩阵提取出来。
out=v$E

可以看到使用limma包的voom函数转化之后，数据更加规整，不会出现数据差距很大的情况。
那整理好数据之后，我们接下来就开始进行cibersort分析，来计算各个样本的免疫细胞含量的情况。
首先我们先把我们这个矩阵导出来，这样也方便我们在后面直接在cibersort函数里面直接读取它。
out=rbind(ID=colnames(out),out)write.table(out,file="symbol.txt",sep="\t",quote=F,col.names=F)

那接下来我们需要去cibersort的官网的下载界面https://cibersort.stanford.edu/download.php去下载22种免疫细胞的注释文件(LM22)以及cibersort的源代码。
我们可以直接在这个界面下载LM22和cibersort的代码，当然小编已经给大家放在附件区，只需要在挑圈联靠公众号回复文末的代号即可获得。
通过网站下载的cibersort代码以及ref.txt的参考文件，还有刚刚我们写出来的表达矩阵，我们接下来就可以进行cibersort函数进行计算每个样本的免疫细胞的含量啦。
首先我们把官网下载cibersort代码使用source函数加载进来。
source("CIBERSORT.R")

接下来按照cibersort函数的要求，分别书写参考文件名和表达矩阵名就可以啦，然后重复1000次。
results=CIBERSORT("ref.txt", "symbol.txt", perm=1000, QN=TRUE)

可以看到经过cibersort计算之后，我们得到了22种免疫细胞含量的表达矩阵，那接下来我们对胃癌病人的22种免疫细胞含量进行柱状图的绘制。
进阶
首先我们先把行名变成第一列，这样才方便我们在后面进行长宽数据转换，同时我们只提取前23列的免疫细胞的数据。
results = as.data.frame(results) %>% rownames_to_column("sample")results = results %>% dplyr::select(1:23)

由于现在的数据是宽数据，我们无法进行ggplot作图，所以我们把宽数据变为长数据，关于长宽数据的转换，先锋宇助教也在之前的推文给大家介绍过三种方法，这里使用我非常喜欢的tidyr包来进行转换。
rt <- pivot_longer(results,-c("sample"), names_to = "immune_cell", values_to = "content")

可以看到我们把刚刚有23列的宽数据变成了3列的长数据，那接下来我们就可以使用ggplot2来进行箱式图的绘制啦。
我们直接把免疫细胞作为X轴，然后免疫细胞含量的content列作为Y轴，然后填充的颜色为免疫细胞那一列。
ggplot(rt, aes(x = immune_cell,y = content,fill = immune_cell)) + geom_boxplot()

可以看到一张绚烂多姿，可以用于发表的箱式图就绘制好啦。
写在最后
先锋宇助教以一张丰富多彩的箱式图提前给大家拜个早年，希望大家在新的一年能够在学习和生活上能够多姿多彩，科研也能做出多姿多彩的成果！！！
2020年小感悟
2020年最大的收获应该是遇见解螺旋，自己也在解螺旋这个平台里面认识了很多优秀的师兄师姐，我从他们身上学习到了很多，无论是对于科研还是未来的人生奋斗目标。最后我想以一句我非常喜欢的曾国藩的一句话来继续迎接充满挑战和希望的2021年！
物来顺应，未来不迎，当时不杂，既过不恋。----《曾国藩家书》
希望大家都能做做好当下的每一件事，不要去忧患未来，对于过去发生的事情就不再留恋。
新春
寄语
希望大家在牛年都能不忘初心，继续为国家的科研事业贡献自己的力量！我们牛年再见，一起牛(扭)转乾坤！
往期传送门：
让生信工作者失业的神器——DrBioRight，真舍不得告诉你！
高效数据清洗！这个R包太强大了！你一定要试试！（文末附赠小彩蛋）
一站式分析R包来了，承包了生信各种分析！太全能了！
搞定这一步，说明你学R有天赋！TCGA数据从下载到差异分析（附代码）
别走，baby，COX森林需要你
生存曲线居然能够批量绘制了
ROC居然能够批量筛选
TCGA转录组和临床信息联合分析，结果发现......
小白重头学习基因组学分析，从绘制棒棒糖图开始
生存曲线一个循环可以画太太太太太太多多多多多多图了吧
华晨宇和张碧晨的瓜吃完了，那基因组和临床信息联合的瓜不来看看吗？
好看的皮囊千篇一律， 高效的R包万里挑一
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
—
END
—

撰文丨先锋宇
排版丨四金兄
值班 | 阿   琛

主编丨小雪球
新年快乐
2020

感谢所有小伙伴的一路陪伴

开心这一路和大家共同成长

2021

我们仍要一起并肩前行

朝更新的目标一起努力


为了感谢大家一路的支持

在春节大年初五迎财神时

酸谈将进行一场
福利抽奖直播
纯抽奖part、全新福利周边
大家一定记得来观看直播奥



直播信息
直播时间：大年初五
直播地点：B站解螺旋直播间

直播内容：福利直播抽奖party

直播地址：
https://live.bilibili.com/8116225
扫码直达直播间



大年初五

不见不散
长按识别二维码免费包邮领取！
继续阅读
阅读原文