找出TCGA中的配对样本并正确展示数据

首先TCGA中的肿瘤样本远大于癌旁样本，甚至有的肿瘤干脆就没有癌旁。

既然有癌旁，那么就有与之对应的癌组织，把他们找出来，可以做很多事情。

开始工作，加载以前的数据，这个数据在我的教程中出现过多次

1. rm(list = ls())

2. load(file = "TCGA_exprSet_plot.Rdata")

3. test <- exprSet[1:10,1:10]

 这个数据是行是样本，行名是TCGA的barcode，列是分组信息和基因名称。

就光这个格式的数据框就可以做很多事情。

甚至牛逼哄哄的协同犯罪(guilt of association)注释基因功能也不在话下

---

下面我们把他们的配对信息找出来

• dplyr 这个包可以对数据框进行增删改查

• tibble 这个包在这里只是方便地实现行名变列，列变行名

• stringr 这个包我们只用到一个功能，就是截取字符串。

配对样本来自于同一个患者，那么TCGA id的前12个字符就是一样的，我们就是根据这个信息来提取配对样本。

1. library(dplyr)

2. library(tibble)

3. library(stringr)

首先提取肿瘤组织的样本,同把截取后的TCGA id作为行名，distinct用于去重，因为有的患者癌症取了多个组织。

1. exprSet_Tumor <- exprSet %>%

2.  rownames_to_column(var="TCGA_id") %>%

3.  filter(sample=="Tumor") %>%

4.  mutate(new=str_sub(TCGA_id,1,12)) %>%

5.  distinct(new,.keep_all = T) %>%

6.  column_to_rownames(var = "new")

看一下获得的数据

1. test <- exprSet_Tumor[1:10,1:10]

可以和一开始数据框比较一下，我把截取后的TCGA id作为行名，为了方便下面提取子集

接下来我们再把癌旁组织信息取出来

1. exprSet_Normal <- exprSet %>%

2.  rownames_to_column(var="TCGA_id") %>%

3.  filter(sample=="Normal") %>%

4.  mutate(new=str_sub(TCGA_id,1,12)) %>%

5.  distinct(new,.keep_all = T) %>%

6.  column_to_rownames(var = "new")

这时候我们需要用intersect获取癌症和癌旁共有的患者,再用这个信息分别调取对应的癌症和癌旁 共有111对配对样本，

1. index <- intersect(rownames(exprSet_Tumor),rownames(exprSet_Normal))

2. exprSet_pair <- data.frame(ID=rep(seq(1,length(index)),2),

3.                           rbind(exprSet_Tumor[index,],exprSet_Normal[index,]))

4. test <- exprSet_pair[,1:10]

第一类，我还增加了配对信息，所以从下到下是1到111，然后又是1到111，总共222行。

---

下面就可以作图展示了，我们选取FOXA1来看一下, R 包就是被很多人追捧的ggpubr。

1. data <- exprSet_pair

2. library(ggpubr)

3. ggpaired(data, x = "sample", y = "FOXA1",

4.         color = "black",

5.         fill = "sample",

6.         line.color = "gray", line.size = 0.4,

7.         ylab = "expression of FOXA1 (VST)",

8.         palette = "npg")+

9.  stat_compare_means(paired = TRUE)

我们看到配对样本是有两条线连接起来的，所以十分方便直观。

除此之外，我们还可以用之前介绍过的方法来作图

1. df <- data[,c("ID","sample","FOXA1")]

2. class <- data$sample

3. pv = wilcox.test(df[,"FOXA1"] ~ class,paired = T)$p.value

4. library(ggplot2)

5. library(cowplot)

6. pv.lab.1 <- "Wilcoxon Signed-Rank Test"

7. pv.lab.2 <- paste("Cancer vs Normal p =", round(pv, 3))

8. lab <- paste("FOXA1\n", pv.lab.1, pv.lab.2, sep="\n")

9. label <- ggdraw() + draw_label(lab)

10. p1 <-ggplot(df, aes(sample, FOXA1,fill=sample)) +

11.  geom_boxplot() +

12.  geom_point(aes(colour=sample), size=2, alpha=0.5) +

13.  geom_line(aes(group=ID), colour="black", linetype="11") +

14.  xlab("") +

15.  ylab(paste("expression of ","FOXA1"," (VST)"))+

16.  theme(legend.position = "none")

17. plot_grid(label,p1,ncol=1, rel_heights=c(.2, 1))

是得，我们可以用别人的包，也可以根据自己的喜好来定制。

现在当我们看这个数据的时候，会发现，这个连接线的方向是不一致的，上面斜着向上，还有一部分样本是斜着向下。 我们肯定会想到，这要看看亚型的信息了。

看亚型我们是有方法的，在这里。

但是，既然我们今天第一次用ggpubr，我们可以开阔一下眼界，我就现学现用

1. ## 亚型

2. # Groups that we want to compare

3. my_comparisons <- list(

4.  c("Basal", "Normal"),

5.  c("Her2", "Normal"),

6.  c("LumA", "Normal"),

7.  c("LumB", "Normal")

8. )

9. # Create the box plot. Change colors by groups: Species

10. # Add jitter points and change the shape by groups

11. ggboxplot(

12.  exprSet, x = "subgroup", y = "FOXA1",

13.  color = "subgroup",

14.  add = "jitter"

15. )+

16.  stat_compare_means(comparisons = my_comparisons, method = "t.test")

这里就一目了然了，确实FOXA1在不同亚型的表现不一样。

我用同一个数据，写了5个以上的帖子，实际上就是阐明一个观点，

尽管漂亮的图表是刚需，无论基金还是文章都需要，却是初学者最不应该学的。初学者应该把有限的精力用在数据清理上。R语言的特点就是R包多，我们只要把数据调整到R包需要的格式，就可以轻易做出漂亮的图片。

我们一开始常常会执着于漂亮的图片，可是最终发现，拥有属于自己的数据才是最重要的。

在这个系列里面，只要我们获取到了正确的数据，然后再用R语言把他调整成一开始的样子，所有剩下的工作就是调包，调包，调包。

稍有能力的朋友，会嘲笑初学者叫调包侠，我反倒是鼓励大家成为调包侠，但是在那之前，最应该掌握的还是数据处理的能力，这在哪里都是基本功。

我如果培训别人也是基于这个思路：