两个CNV信息的合并后比较

答读者问：

有朋友问到他对同一个样本SNP6.0芯片测到了CNV信息，也做了WES得到了CNV信息，该如何比较这两个CNV结果呢？（PS：其实他提问非常含糊，描述也很差，但考虑到对方是生信技能树的VIP，我还是硬着头皮把这个问题提炼出成为了一个可解决的。）

很多情况下，CNV-calling工具对不同测序数据的bam文件得到的segment信息坐标是基本上完全不一致的。

如下所示：

  sample   chr     start       end       cnv
1 MR1.bam chr10   3259559 130018550 -0.508188
4 MR1.bam chr11   3131777  57645966 -0.525031
5 MR1.bam chr11  57721045  69579301  0.402708
7 MR1.bam chr11  69591008 101526888 -0.416663 

也就是说不同样本的cnv片段的染色体起始和终止坐标肯定是没办法完全一样，而且不同的工具或者测序手段得到的分辨率不一样，也很难直接比较。

这里有两种方法来处理，首先可以以基因或者外显子来进行映射，这样可以统一坐标，但是bin太小了，而且忽略了基因组大部分区域。

另外一种方法是需要分区间来统一，比如1MB的区间：

head(a) 
# 上面的文件内容存储在a这个变量里面 
colnames(a)=c('sample','chr','start','end','cnv')
## 因为以1Mb为区间，就简单粗暴的删除小片段CNV
a=a[a$end-a$start > 1000000,]
tmp = apply(a, 1, function(x){
  #x=as.character(a[1,])
  tmp = lapply(seq(round(as.numeric(x[3])/1000000),
             round(as.numeric(x[4])/1000000),1), function(i){
     return(c(x[1],paste(x[2],i,sep=":"),x[5]))
  })
  tmp=do.call(rbind,tmp)
  return(tmp)
})
cnv_wes=do.call(rbind,tmp)
head(cnv_wes)

转换后的外显子的CNV信息如下：

    sample      pos       cnv        
[1,] "MR1" "chr10:3" "-0.508188"
[2,] "MR1" "chr10:4" "-0.508188"
[3,] "MR1" "chr10:5" "-0.508188"

可以看到之前的chr,start,end 被 统一为了pos列，值得注意的是原来的拷贝数变异的segment文件每个样本只有几百行，归一化后就是染色体的长度啦，应该是有30G/1MB, 就是每个样本应该是3000行 ！

更加值得注意的是，上面的代码并非是完美的，希望大家能交流自己的方法，谢谢。

总之，有了统一的坐标，这样就可以合并做热图咯！

D[1:10,1:10]
pos$pos=paste0("chr",pos$chr,":",round(pos$start/1000000))
### 对那些比 1Mb小的CNV区间，就需要另行处理，如下：
d_n=apply(D,1,function(x){
  #x=as.numeric(D[2,])
  tmp = aggregate(x, list(pos=pos$pos) ,mean)
  return(tmp[,2])
})
## 这里是取平均值。
x=as.numeric(D[1,])
tmp = aggregate(x, list(pos=pos$pos) ,mean)
rownames(d_n)=tmp[,1]

如果要把两个矩阵合并，通常是需要 保证大家的pos是一致的，示例代码如下：

table(rownames(d_n) %in%   cnv_wes_pos )
table( cnv_wes_pos  %in%   rownames(d_n) )
pos=intersect( cnv_wes_pos,rownames(d_n))

library(stringr)
pos_df=str_split(pos,':',simplify = T)
pos_df=as.data.frame(pos_df)
colnames(pos_df)=c('chr','start')
pos_df[,2]=as.numeric(pos_df[,2])
pos_df$chr=paste('chr',
                 sprintf('%02d',as.numeric(gsub('chr','',pos_df$chr))) ,
                 sep = '') 

od=order(pos_df[,1],pos_df[,2])
pos_df=pos_df[od,]
pos=pos[od]


head(pos_df)
head(pos)


table(cnv_wes[,1])

上面的代码只是一个示例，具体情况需要根据自己的理解来修改。

而且应该是有更好的方法来解决，比如分辨率过高的CNV信息也可以使用segment算法先片段化，而不是简单粗暴的在1Mb区间内取平均值。

欢迎大家提出自己的想法，或者email 来信与我交流。