项目实战 |根据找到的variants的结果生成突变矩阵

最新推荐文章于 2022-11-06 10:01:05 发布

今天也是个妖精头子呀

最新推荐文章于 2022-11-06 10:01:05 发布

阅读量772

点赞数

分类专栏：谱系追踪

本文链接：https://blog.csdn.net/weixin_40640700/article/details/118636507

版权

谱系追踪专栏收录该内容

82 篇文章 12 订阅

订阅专栏

这个代码要解决的具体需求如下：
（1）将每一个细胞的突变矩阵读入R中（主要有46个数据矩阵）
（2）筛选突变位点（要求：至少在4个细胞中出现）
（3）整合生成最终的突变位点集
（4）计算每一个突变位点在每一个细胞中对应的VAF值
（5）整合各个细胞的VAF值，形成突变矩阵
（6）根据突变矩阵，绘制heatmap图，完成任务

我们需要我们的程序能够完成这些事儿，现在开始（16：06）。

参考链接：https://cloud.tencent.com/developer/article/1763129
发现了一个叫tidyverse的包，使用这个包，我们可以很方便的读取这些数据信息。

setwd("F://Rworkplace//mutation") #数据集所在的位置
files<-dir(path = "F://Rworkplace//mutation",full.names = T,pattern = ".txt") #读取目录下的所有数据
library(tidyverse) #加载包
df<-map(files,read.table,sep="\t") #批量的读取数据文件
class(df) #这个包的功能就是把数据读取到数据集中

“list”

离支线太远了（17：58）。我总是很容易，钻牛角尖。现在同一个问题，弄了两个小时了还没有推进一点点，心态真的很糟糕。
越着急越不行，毕竟这个问题的确超出了我的能力，唯有静下心来，心平气和的，才有突破的可能。

目前试过两种方法：
（1）reduce（intersect，list）

报错为：Error in x[[1]] : object of type ‘closure’ is not subsettable
而换一个list则完全可以使用这个函数，我想不明白为什么？

（2）upsetR

报错为：1：ncol() zero
我想不明白这个函数的使用方法，明明提交的也是一个list，为什么运行不了？

所以，综合使用这两种方法，都不能够解决我现在的问题，我现在该怎么办？
是再去百度上检索，还是自己静下心来想一想，如果按照自己目前掌握的方法，我要怎么处理这件事？

我可以这个样子处理：
（1）首先把这些位点全部的取出来，不去重的那种，

#取出有start坐标的那一行
sublist<-function(x){
  as.numeric(unlist(x$V2[-1]))
}
listInput<-sapply(df, sublist)
length(listInput)
#我感觉对于我这种小菜鸟，用循环来批量处理是最现实的
b<-c()
for(i in (1:length(listInput))){
  a<c()
  a<-as.numeric(unlist(listInput[[i]]))
  b<-append(b,a)
}
site_freq<-as.data.frame(table(b))
#这个list就是我们要的，我们可以根据这个list来筛选位点
dim(site_freq)
[1] 67999     2

得到的site_freq的数据矩阵：
在这里插入图片描述
（2）然后我们可以根据这个频次矩阵，从中筛选出我们想要的位点。
筛选标准：

frequency ≥ 3
frequency=46的位点删去（也就是说在所有的文件中都突变，没有比较的意义）

filter_site<-subset(site_freq,site_freq$Freq>=3&site_freq$Freq!=46)
dim(filter_site)
[1] 7488    2

（3）然后将这些位点，存入到一个vector中。
我们来看一下这个数据矩阵。

# interest_site<-as.numeric(filter_site[,1]) #这样写是错误的
interest_site<-filter_site[,1]
#我终于明白哪里错了,我想要取b列的值，但是呢！它取成了行号
#所以为什么？
#应该是原先的as.numeric()出了问题
#length(intersect(interest_site,df[[1]]$V2))
#mutation_cell_1<-subset(df[[1]],subset=(df[[1]]$V2%in%interest_site))
#我的想法很简单，就是生成多个变量就行，一步一步来

（4）从原始的各个数据集中，取出start位点属于这个数据集的位点的行。
这里需要用到一个属于的判断符号%in%。我们来进行处理。

for (i in 1:length(df)) {
  assign(paste("mutation_cell",i,sep = "_"),subset(df[[i]],subset=(df[[i]]$V2%in%interest_site)))
} #这里的assign，我研究了好久，可以作为经验写出来 #包括上次拆分那个数据是同样的道理

#这个循环的嵌套关系比较复杂，但是从本质上来讲，就是说将df中的每一个list的值都取出来，然后逐个的进行判断。最后得到结果data.frame尝试赋值给循环定义的object。

在这里插入图片描述
可以大致的看到这些细胞中最多含有的符合要求的variant的数目不过是十几个。我觉得这个筛选的过程是出错的。
因为我们前面满足要求的filter_site有7000左右，并且呢！我看了一下这个variant中好像只有chrM的结果，肯定是有问题的。
我们再次检查一下我们的代码。

经过检查代码，发现，是我们在取感兴趣的位点的时候，程序错误理解为取行号，所以导致最后能够匹配的数量很少。所以，回到上面，进行了调整。

调整后继续匹配，查看最后的结果。
在这里插入图片描述
这样明显多多了！是契合我们的预期的。我们打开一个数据矩阵看一下！找到了常染色体上的突变位点信息。

所以，到这一步骤，我们check的结果都是正确的。

（5）随后，我们计算每一个位点的VAF值，对于每一个细胞，最后产生两列的data.frame，第一列是由原始的前面五列就是位点的信息缩合的，最后一列是我们计算得到的这个细胞在这个位点的VAF值。

今天或许可以把这个层次的突变矩阵画出来。关于observable的位点，我们稍后再进行筛选。
突然有些感慨，问题的解决方案的发现，有时候是需要一些运气的因素的。有时候，纠缠于某个问题挑不出来，一个偶然的注意，撬开了一点点的思路，继续往前挖。也许就是通天大道。
我还是非常欣赏胡适先生的那句话的，大胆假设，小心求证。不要急急忙忙的寻找问题的结果，而是一点点的去证实它的可靠性，这一步确保了，才能说下一步是有意义的。
我相信白老师答应我了，就一定会践行承诺的。我想，所以，我不必着急，仅仅是实事求是的汇报就好了。完成比完美更重要，先把事情完成，再去追求完美。
我觉得认清自己很重要，我就是一个啥基础都没有的实习生，是努力的在做我应该做的事情了，做不好是很正常的事情，承认自己比较菜，会有一种坦然。
但是在菜的前提下，不要忘记认真努力才是啊！

在做之前，首先明确VAF值应该如何去求？
这里其实也遇到了一个问题（就是我昨天一直试图想要去解决的）：

我们使用的是tumor_only的方式（也就是说是没有对照的normal样本的），这样calling出来的结果在覆盖度方面比较奇怪。
我前前后后有查看mtDNA和我现在用CML样本跑的样本，我发现用tumor-only跑的样本的coverage的参数和我们在vcf文件中使用ref+alt的值不太一样，有时候甚至会发现在coverage中位点为零的位点会出现在vaf的样本中，就让我产生了比较大的疑惑。
我感觉做我们这一行，一定要非常精确才是，我们要对我们自己的结果负责，只有我们能力强了，我们的结果才有可能更加的精确。我感觉科学就是一个不断的怀疑与自我怀疑的过程。

我现在先计算VAF值，将其缩写成一个两列的data.frame。

计算VAF值，并产生VAF列。

for (i in 1:dim(mutation_cell_1)[1]){
  mutation_cell_1$VAF[i]<-as.numeric(mutation_cell_1$V10[i])/(as.numeric(mutation_cell_1$V9[i])+as.numeric(mutation_cell_1$V10[i]))
  }

注意将其as.numeric转换，转换之后的可以进行计算，否则会出错。

Error in mutation_cell_1 $V9[1] + mutation_cell_1$ V10[1] :
non-numeric argument to binary operator

现在计算完成。

然后将variant的label转换为“–>”

mutation_cell_1<-within(mutation_cell_1, mutation<- paste(V4,V5,sep="->"))

对数据进行提取，合并

mutation_cell_1_VAF<-cbind(mutation_cell_1$V1,mutation_cell_1$V2,mutation_cell_1$mutation,mutation_cell_1$VAF)

（6）现在尝试将上面这一步，扩展到这46个数据中，生成46个VAF矩阵。
生成函数。

cal_VAF<-function(data,m){
  for (i in 1:dim(data)[1]){
    data$VAF[i]<-as.numeric(data$V10[i])/(as.numeric(data$V9[i])+as.numeric(data$V10[i]))
  }
  data<-within(data, mutation<- paste(V4,V5,sep="->"))
  assign(m,cbind(data$V1,data$V2,data$mutation,data$VAF))
  }

46个数据文件生成（虽然使用了函数，可是输入变量的时候依旧是很麻烦）。

mutation_cell_2_VAF<-cal_VAF(mutation_cell_2,"mutation_cell_2_VAF")
mutation_cell_3_VAF<-cal_VAF(mutation_cell_3,"mutation_cell_3_VAF")
mutation_cell_4_VAF<-cal_VAF(mutation_cell_4,"mutation_cell_4_VAF")
mutation_cell_5_VAF<-cal_VAF(mutation_cell_5,"mutation_cell_5_VAF")
mutation_cell_6_VAF<-cal_VAF(mutation_cell_6,"mutation_cell_6_VAF")
mutation_cell_7_VAF<-cal_VAF(mutation_cell_7,"mutation_cell_7_VAF")
mutation_cell_8_VAF<-cal_VAF(mutation_cell_8,"mutation_cell_8_VAF")
mutation_cell_9_VAF<-cal_VAF(mutation_cell_9,"mutation_cell_9_VAF")
mutation_cell_10_VAF<-cal_VAF(mutation_cell_10,"mutation_cell_10_VAF")
mutation_cell_11_VAF<-cal_VAF(mutation_cell_11,"mutation_cell_11_VAF")
mutation_cell_12_VAF<-cal_VAF(mutation_cell_12,"mutation_cell_12_VAF")
mutation_cell_13_VAF<-cal_VAF(mutation_cell_13,"mutation_cell_13_VAF")
mutation_cell_14_VAF<-cal_VAF(mutation_cell_14,"mutation_cell_14_VAF")
mutation_cell_15_VAF<-cal_VAF(mutation_cell_15,"mutation_cell_15_VAF")
mutation_cell_16_VAF<-cal_VAF(mutation_cell_16,"mutation_cell_16_VAF")
mutation_cell_17_VAF<-cal_VAF(mutation_cell_17,"mutation_cell_17_VAF")
mutation_cell_18_VAF<-cal_VAF(mutation_cell_18,"mutation_cell_18_VAF")
mutation_cell_19_VAF<-cal_VAF(mutation_cell_19,"mutation_cell_19_VAF")
mutation_cell_20_VAF<-cal_VAF(mutation_cell_20,"mutation_cell_20_VAF")

##########################################################################

mutation_cell_21_VAF<-cal_VAF(mutation_cell_21,"mutation_cell_21_VAF")
mutation_cell_22_VAF<-cal_VAF(mutation_cell_22,"mutation_cell_22_VAF")
mutation_cell_23_VAF<-cal_VAF(mutation_cell_23,"mutation_cell_23_VAF")
mutation_cell_24_VAF<-cal_VAF(mutation_cell_24,"mutation_cell_24_VAF")
mutation_cell_25_VAF<-cal_VAF(mutation_cell_25,"mutation_cell_25_VAF")
mutation_cell_26_VAF<-cal_VAF(mutation_cell_26,"mutation_cell_26_VAF")
mutation_cell_27_VAF<-cal_VAF(mutation_cell_27,"mutation_cell_27_VAF")
mutation_cell_28_VAF<-cal_VAF(mutation_cell_28,"mutation_cell_28_VAF")
mutation_cell_29_VAF<-cal_VAF(mutation_cell_29,"mutation_cell_29_VAF")


mutation_cell_30_VAF<-cal_VAF(mutation_cell_30,"mutation_cell_30_VAF")
mutation_cell_31_VAF<-cal_VAF(mutation_cell_31,"mutation_cell_31_VAF")
mutation_cell_32_VAF<-cal_VAF(mutation_cell_32,"mutation_cell_32_VAF")
mutation_cell_33_VAF<-cal_VAF(mutation_cell_33,"mutation_cell_33_VAF")
mutation_cell_34_VAF<-cal_VAF(mutation_cell_34,"mutation_cell_34_VAF")
mutation_cell_35_VAF<-cal_VAF(mutation_cell_35,"mutation_cell_35_VAF")
mutation_cell_36_VAF<-cal_VAF(mutation_cell_36,"mutation_cell_36_VAF")
mutation_cell_37_VAF<-cal_VAF(mutation_cell_37,"mutation_cell_37_VAF")
mutation_cell_38_VAF<-cal_VAF(mutation_cell_38,"mutation_cell_38_VAF")
mutation_cell_39_VAF<-cal_VAF(mutation_cell_39,"mutation_cell_39_VAF")
mutation_cell_40_VAF<-cal_VAF(mutation_cell_40,"mutation_cell_40_VAF")


mutation_cell_41_VAF<-cal_VAF(mutation_cell_41,"mutation_cell_41_VAF")
mutation_cell_42_VAF<-cal_VAF(mutation_cell_42,"mutation_cell_42_VAF")
mutation_cell_43_VAF<-cal_VAF(mutation_cell_43,"mutation_cell_43_VAF")
mutation_cell_44_VAF<-cal_VAF(mutation_cell_44,"mutation_cell_44_VAF")
mutation_cell_45_VAF<-cal_VAF(mutation_cell_45,"mutation_cell_45_VAF")
mutation_cell_46_VAF<-cal_VAF(mutation_cell_46,"mutation_cell_46_VAF")

明天，将这46个文件，整合到一起。

唐老师的数据先放一放，我先尝试着接着处理完接下来的数据，争取在今晚之前（18:00）绘制好热图，明天阶段性的向老师汇报。先做最紧急最重要的事情。我觉得当老师还挺麻烦的，需要多线程的处理那么多的事情。
现在开始。（14:10）

（7）多个样本，根据公共的行，整合在一起。

我们本次使用的还是reduce（）这个函数，总的来说，虽然样子丑了点，但是可以解决我们的问题。

filter_site_update<-as.data.frame(filter_site[,-2])
colnames(filter_site_update)<-"V2"
#merge(filter_site_update,mutation_cell_1_VAF,by.x= "V2",by.y = "V2",all = TRUE)
#再次尝试使用reduce函数
test<-Reduce(function(x, y) merge(x, y, all=TRUE,by.x="V2",by.y="V2"),
             list(filter_site_update,
                  mutation_cell_1_VAF,mutation_cell_2_VAF,mutation_cell_3_VAF,
                  mutation_cell_4_VAF,mutation_cell_5_VAF,mutation_cell_6_VAF,
                  mutation_cell_7_VAF,mutation_cell_8_VAF,mutation_cell_9_VAF,
                  mutation_cell_10_VAF,mutation_cell_11_VAF,mutation_cell_12_VAF,
                  mutation_cell_13_VAF,mutation_cell_14_VAF,mutation_cell_15_VAF,
                  mutation_cell_16_VAF,mutation_cell_17_VAF,mutation_cell_18_VAF,
                  mutation_cell_19_VAF,mutation_cell_20_VAF,mutation_cell_21_VAF,
                  mutation_cell_22_VAF,mutation_cell_23_VAF,mutation_cell_24_VAF,
                  mutation_cell_25_VAF,mutation_cell_26_VAF,mutation_cell_27_VAF,
                  mutation_cell_28_VAF,mutation_cell_29_VAF,mutation_cell_30_VAF,
                  mutation_cell_31_VAF,mutation_cell_32_VAF,mutation_cell_33_VAF,
                  mutation_cell_34_VAF,mutation_cell_35_VAF,mutation_cell_36_VAF,
                  mutation_cell_37_VAF,mutation_cell_38_VAF,mutation_cell_39_VAF,
                  mutation_cell_40_VAF,mutation_cell_41_VAF,mutation_cell_42_VAF,
                  mutation_cell_43_VAF,mutation_cell_44_VAF,mutation_cell_45_VAF,
                  mutation_cell_46_VAF))

得到的数据矩阵比较丑。如下图。

在这里插入图片描述
（8）上面得到的这个数据矩阵，距离我们的“理想”矩阵还是有些区别的。所以，接下来就是对于这个数据矩阵进行“再整理”。

将行名转化为类如chr1 100292207 T->C这种形式。
每一列的列名，变为cell1、cell2、cell3、cell4等等。
整理成，行是每一个突变位点的突变类型，列为每一个细胞，中间的每一个值是一个关于每一个细胞在每一个突变位点的VAF值的这样的突变矩阵。

dplyr包有时间可以多学一下。这个包根据这几天的了解看，是用来整理数据框的，可以说是data.frame的进阶版，我自己倒挺有兴趣的。

在这个过程中遇到了一个问题，对这个数据的dm进行检测的时候，发现我们生成的这样的一个数据矩阵的行数大于我们的位点数据。追踪其原因是由于，对于某些特定的细胞，其突变位点对应于两种类型的突变。

dup<-variant_data[duplicated(variant_data[,1]),1]
#FALSE  TRUE 
#7488   531

这个问题先暂时放一下，我现在的想法是先解决这个行名的问题，把关于突变的label的列聚在一起。
尝试对染色体号进行去重的过程中，又出现了新的问题。

和我原来预期的不太一样的是，我以为在染色体上，某一个位点是确定且唯一的。但是，这个结果就让人很惊讶。而且，相应位点的对应的碱基也不同。
chrM 2805 A 48
chrUn_KI270582v1 2805 T 1
chrUn_GL000218v1 2805 A 1
chrUn_KN707706v1_decoy 2805 a 1
chrUn_KN707896v1_decoy 2805 c 13

那遇到这种情况该怎么办？该怎么理解这种现象呢？
同样是忽略这种情况。

f<-c()
for (j in 1:dim(variant_data_1)[1]){
  d<-c()
  for (i in 2:47){
    if(!is.na(variant_data_1[j,i])){
      b<-c()
      b<-variant_data_1[j,i]
      d<-append(d,b)
  }else{
      next
   # print(d)
   # print(j)
  }
    #len<-length(unique(d))
    #f<-append(f,len)
    if(length(unique(d))==1){
      variant_data_1[j,2]<-unique(d) #问题出在这里
    }
  }
}

就是在这一步将这些位点剔除（肯定是在数据筛选的过程中会有一些损耗的）。
下一步就是将其他的第3~46列剔除，保留第2列。

variant_data_2<-variant_data_1[,-c(3:47)]

剔除之后，继续按照同样的方法，选取突变类型。

f<-c()
for (j in 1:dim(variant_data_2)[1]){
  d<-c()
  for (i in 3:48){
    if(!is.na(variant_data_2[j,i])){
      b<-c()
      b<-variant_data_2[j,i]
      d<-append(d,b)
    }else{
      next
      # print(d)
      # print(j)
    }
    #len<-length(unique(d))
    #f<-append(f,len)
    if(length(unique(d))==1){
      variant_data_2[j,3]<-unique(d) #问题出在这里
    }
  }
}
variant_data_3<-variant_data_2[,-c(4:48)]

这个事情，处理结束之后（肯定也会剔除一些东西，但是剔除多少我没有统计）。
紧接着，我们将几列合并，生成我们想要的突变类型的列。这个我们之前有处理过，我来看一下。

ariant_data_4<-within(variant_data_3, mutation_site<- paste(V1.x,V2,V3.x,sep="\t"))
dim(variant_data_4)
variant_data_5<-variant_data_4[,c(50,4:49)]

接着，对行和列进行重命名。
首先生成一个字符串的序列。

cell_label<-paste("cell",1:46,sep = "_")
colnames(variant_data_5)<-c("variant",cell_label)
#果然，如我猜测的那样，有重复的mutation
#出现这些冗余行的根本原因，要么是位点的重复，要么是同一个突变位点，有两个calling结果
rownames(variant_data_5)<-variant_data_5[,1]

这个时候出现了一些问题：

Error in .rowNamesDF<-(x, value = value) : 不允许有重复的’row.names’
In addition: Warning message:
non-unique values when setting ‘row.names’: ‘chr1 40072268 A->C’, ‘chr10 3776932 C->CT’, ‘chr11 27496002 C->CA’, ‘chr13 52699827 G->GT’, ‘chr14 59371109 C->A’, ‘chr16 66565991 GA->G’, ‘chr17 44223561 A->AC’, ‘chr17 78225123 G->A’, ‘chr20 49094159 G->A’, ‘chr3 177021370 CA->C’, ‘chr4 41960765 G->GA’, ‘chr5 134750980 C->CA’, ‘chr6 143830177 T->TA’, ‘chr7 26212582 T->G’, ‘chr8 25509718 TAA->T’, ‘chrX 25645675 T->A’

所以，如我之前的料想，出现了一些冗余的行。
主要原因，要么是位点的重复，要么是突变类型的重复。
现在的想法就是全部的剔除这些行（甚至一个都不保留，因为这些是在处理的过程中ambiguous的部分）

ambiguious_site<-duplicated(variant_data_5[,1])
ambiguious_variant<-variant_data_5[ambiguious_site,1]

variant_data_6<-variant_data_5[!ambiguious_site,] #剔除了ambiguous的位点（但仍然保留唯一的位点，这种方式私以为不够准确）

rownames(variant_data_6)<-variant_data_6[,1]
variant_data_7<-variant_data_6[,-1]
rownames(variant_data_7)<-variant_data_6[,1] #对数据进行重新的整理

write.table(variant_data_7,"variant.txt",quote = F)

#现在先做一个处理，NA用0来表示
variant_data_7[is.na(variant_data_7)] <- 0
variant_data_7<-apply(variant_data_7,2,as.numeric)
library(pheatmap)
pheatmap(variant_data_7,cluster_rows = T,cluster_cols = F,color = colorRampPalette(c("red","yellow","white"))(100),show_colnames = F) #绘制heatmap

variant_chr6_A<-variant_data_6["chr6	132814747	A->G",-1]
variant_chr6_A[is.na(variant_chr6_A)]<-0

breaks<-seq(0,1,length.out =20) #之所以出错是因为break在R中是有意义的，不能用于变量名
hist(as.numeric(variant_chr6_A),breaks= c(0.00,0.125,0.25,0.375,0.50,0.625,0.75,0.875,1.00),ylim = c(0,50),xlim =c(0,1),xlab = "Variant Allele Frequence",ylab="Count",main = "chr6 \t 132814747	A->G")


variant_chr6_G<-variant_data_6["chr6	35469489	G->GA",-1]
variant_chr6_G[is.na(variant_chr6_G)]<-0

breaks<-seq(0,1,length.out =20) #之所以出错是因为break在R中是有意义的，不能用于变量名
hist(as.numeric(variant_chr6_G),breaks= c(0.00,0.125,0.25,0.375,0.50,0.625,0.75,0.875,1.00),ylim = c(0,50),xlim =c(0,1),xlab = "Variant Allele Frequence",ylab="Count",main = "chr6 \t 35469489	G->GA")

最后生成的结果图如下：
在这里插入图片描述

在这里插入图片描述

完成。

今天也是个妖精头子呀

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
项目实战 |根据找到的variants的结果生成突变矩阵

这个代码要解决的具体需求如下：（1）将每一个细胞的突变矩阵读入R中（主要有46个数据矩阵）（2）筛选突变位点（要求：至少在4个细胞中出现）（3）整合生成最终的突变位点集（4）计算每一个突变位点在每一个细胞中对应的VAF值（5）整合各个细胞的VAF值，形成突变矩阵（6）根据突变矩阵，绘制heatmap图，完成任务我们需要我们的程序能够完成这些事儿，现在开始（16：06）。...
复制链接

扫一扫

专栏目录