两个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区间内取平均值。

两个CNV信息的合并后比较

相关推荐