CNTools 取得 CNV 矩阵
2020-07-03 本文已影响0人
BeeBee生信
CNTools 主要功能是将多个样品 CNV 结果转换成矩阵。
CNTools原理功能
输入数据为 R 包 DNAcopy 的 segment
函数结果。在 TCGA 数据集从 GDC 下载 CNV 结果是 DNAcopy 已经处理好的 seg 格式文件,不需要用 segment
函数处理。这时需要手动将所有的 seg 文件合并,直接 bind_rows
就行,最好把第一列改成 TCGA Barcode 方便后续分析。GDC 下载的 seg 文件一般来说一个样本有 4 文件,2个 Normal 组织2个 Tumor 组织,每一种组织里1个原始 seg 文件1个 Masked seg。后者是在前者基础上进行了过滤,将 Y 染色体和正常组织常见 CNV 移除,所以可以近似认为是 Somatic CNV(Tumor组织时)。 Masked seg 文件名带有 nocnv
几个字,差不多 *.nocnv_grch38.seg.v2.txt
这种格式。
下面是步骤,合并好的 seg 文件要先修改表头。
> colnames(segData) <- c("ID", "chrom", "loc.start", "loc.end", "num.mark", "seg.mean")
> head(segData, 3)
ID chrom loc.start loc.end num.mark seg.mean
1 TCGA-2F-A9KO-01A 1 3301765 28398187 13436 -0.2143
2 TCGA-2F-A9KO-01A 1 28399444 40918052 6652 0.0221
3 TCGA-2F-A9KO-01A 1 40918166 46819458 3203 -0.1991
还需要基因信息表,这个表头和每列内容也是锁死的,按照下面这个案例来。其中 geneid
是 entrez ID。这个表格需要自己去整理出来。
> head(geneInfo)
chrom start end geneid genename
1 19 58345178 58353492 1 A1BG
2 12 9067664 9116229 2 A2M
3 8 18170477 18223689 9 NAT1
4 8 18391282 18401218 10 NAT2
5 14 94612384 94624055 12 SERPINA3
6 3 151814073 151828488 13 AADAC
然后是读取 seg 数据并转换到矩阵的过程。
> cnSeg <- CNSeg(segData)
> rdGene <- getRS(cnSeg, by="gene", imput=FALSE, XY=FALSE, geneMap=geneInfo, what = "mean")
> rsGene <- rs(rdGene)
> head(colnames(rsGene), 8)
[1] "chrom" "start" "end" "geneid"
[5] "genename" "TCGA-2F-A9KO-01A" "TCGA-2F-A9KP-01A" "TCGA-2F-A9KQ-01A"
得到矩阵前几列是基因信息,后面是每个样本该基因的 CNV 情况,在我这例子用 what="mean"
所以是 log2 segment 平均值。