基因组学

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 平均值。

上一篇下一篇

猜你喜欢

热点阅读