内容

方法的所有可用函数，本文档将为您提供指导iGC包中。包iGC旨在同时分析基因表达(GE)和拷贝数改变(CNA)。

传统的CNA分析方法是研究不同类型的样本，并通过维恩图将其结果进行整合。然而，当在多个平台上观察到低再现性和不一致性时，挑战就出现了。为了解决这些问题，iGC同时测试基因表达谱和拷贝数变异。

有关iGC使用的方法的更多信息，请参阅我们的出版物:赖一品，王良博，赖良川，蔡蒙勋，吕子品，Eric Y Chuang。igc -基因表达与拷贝数改变综合分析包Bioinfomatics(即将出版)。

0．1安装

iGC安装在Bioconductor上，可以按照标准安装程序进行安装。

如果(!requireNamespace("BiocManager"， quiet =TRUE) install.packages("BiocManager") BiocManager::install("iGC")

使用,

图书馆(iGC)

0.2一般工作流程

一般工作流程可以总结如下:

基本上有四个步骤，对应四个R函数，来完成分析:

  1. 指定样本及其数据关系
  2. 阅读和组织基因表达文件
  3. 根据基因区域读取和组织CNA文件
  4. 鉴定na驱动的差异表达基因

0.3数据源

如图所示的工作流程，样本配对分析需要CNA和基因表达数据。由于iGC在基因水平或染色体位置读取CNA和基因表达数据，因此数据源是平台和技术不可知的，即来自微阵列或下一代测序的数据都是可接受的。然而，这也意味着需要一些标准的预处理步骤来将原始读数转换为可解释的数据。

0.3.1基因表达

对于每个样本，数据应该至少有两列，基因和表达:

基因表达A 0.1 B -0.5 C 0.4

iGC读取所有样本的基因表达create_gene_exp并制作一个联合基因表达表，

GENE SampleB…A 0.1 0.5 b -0.5 2.1 c 0.4 na

这意味着，如果存在这样的联合表达式表，可以跳过这一步。

0.3.2中央社

iGC接受两种CNA数据格式:基于染色体位置的和基于基因的create_gene_cna而且direct_gene_cna分别。对于基于染色体定位的CNA记录，每个样本的数据至少应该有以下列:

3.染色体开始结束表达式1 1000 5030 2.5 10 10 2560 -1.3 X 12345 14200

iGC将通过查找基因组参考将这种格式的记录转换为基于基因的格式，目前仅支持hg19(参见问答了解更多信息和未来的发展计划)。

对于基于基因的CNA记录，数据格式类似于基因表达数据，

基因表达A 2.5 B 2.5 C -1.3

最后，iGC读取所有样本的CNA记录，生成联合CNA表。与基因表达处理不同，CNA记录根据给定的阈值额外转换为CNA增益/损失事件。

GENE SampleB…A 0 0 b -1 -1 c 0 1

如果这样的联合表已经存在，也可以跳过这一步。

0.3.3自定义阅读器功能

支持各种数据格式，所有create_gene_exp，create_gene_cna,direct_gene_cna接受一个论点read_fun，一个读取示例数据的函数，以自定义如何读取数据。请参阅它们的文档以了解实施细节。

0.4使用的例子

这里我们展示了一个使用iGC附带的数据的示例。

0.4.1数据来源示例

的用法iGC该包装包含50个乳腺癌样本，这些样本是从TCGA 3级总共523个乳腺癌样本的微阵列数据集中选择的。

对于每个样本，它包含一个配对基因表达(GE)和拷贝数(CN)数据。GE数据由Agilent G4502A平台进行，CN数据由Genome Wide SNP平台生成。

0.4.2样本描述生成

首先，创建一个示例描述表，将示例名称与其数据连接起来。它可以存储为CSV文件，有三列:样本，CNA_filepath,GE_filepath．

Sample_desc_pth <- system。file("extdata"， "sample_desc.csv"， package = "iGC")

或者，可以传递三个不同的字符向量来创建同一个表，

sample_desc <- create_sample_desc(sample_names = sample_desc$Sample, CNA_filepath = sample_desc$CNA_filepath, GE_filepath = sample_desc$GE_filepath)
头(sample_desc)
# # # # 1:样本TCGA-AO-A0JL-01A # # 2: TCGA-BH-A0HF-01A # # 3: TCGA-BH-A0HK-01A # # 4: TCGA-AO-A0JF-01A # # 5: TCGA-BH-A0DP-01A # # 6: TCGA-AO-A0JJ-01A # # CNA_filepath # # 1: / tmp / RtmpU2E9r7 / Rinst1c75e61cdf452e / iGC / extdata / CNA A04_697150.hg19.txt # # 2: / tmp / RtmpU2E9r7 / Rinst1c75e61cdf452e iGC / extdata / CNA / A05_697032.hg19.txt # # 3: / tmp / RtmpU2E9r7 / Rinst1c75e61cdf452e iGC / extdata / CNA / A07_697200.hg19.txt # # 4: / tmp / RtmpU2E9r7 / Rinst1c75e61cdf452e iGC / extdata / CNA / A09_697146.hg19.txt # # 5:/ tmp / RtmpU2E9r7 / Rinst1c75e61cdf452e iGC / extdata / CNA / B01_697160.hg19.txt # # 6: / tmp / RtmpU2E9r7 / Rinst1c75e61cdf452e / iGC / extdata / CNA B03_697186.hg19.txt # # GE_filepath # # 1: / tmp / RtmpU2E9r7 / Rinst1c75e61cdf452e iGC / extdata /电气/ US82800149_251976012210_S01.tcga_level3.txt # # 2: / tmp / RtmpU2E9r7 / Rinst1c75e61cdf452e iGC / extdata /电气/ US82800149_251976012504_S01.tcga_level3.txt # # 3: / tmp / RtmpU2E9r7 / Rinst1c75e61cdf452e iGC / extdata /电气/ US82800149_251976012510_S01.tcga_level3.txt # # 4:/tmp/RtmpU2E9r7/Rinst1c75e61cdf452e/iGC/extdata/GE/US82800149_251976012517_S01.tcga_level3.txt ## 5: /tmp/RtmpU2E9r7/Rinst1c75e61cdf452e/iGC/extdata/GE/US82800149_251976012261_S01.tcga_level3.txt ## 6: /tmp/RtmpU2E9r7/Rinst1c75e61cdf452e/iGC/extdata/GE/US82800149_251976012204_S01.tcga_level3.txt ##

0.4.3联合基因表达表

其次，我们将所有样本的基因表达连接成一张表。

gene_exp <- create_gene_exp(sample_desc, progress = FALSE)

create_gene_exp内置读取功能，读取基因表达数据。但是对于无法识别的格式，可以通过定义自己的读取器函数read_fun，

Gene_exp <- create_gene_exp(sample_desc, read_fun = read。table, progress = TRUE, progress_width = 60， #参数传递给自定义read_fun(这里是read.table)头= FALSE, skip = 2, na。string = "null"， colClasses = c("character"， "double"))

注意参数头，跳过，na.srings,colClasses不被使用create_gene_exp但是传给了read.table，这里定义的自定义阅读器函数。

我们选取前9个样本中TP53、BRCA1和NFKB1基因的表达(第一列包含基因名)。

gene_exp[基因%in% c('TP53'， 'BRCA1'， 'NFKB1')， 1:10, with=FALSE]
# #基因TCGA-AO-A0JL-01A TCGA-BH-A0HF-01A TCGA-BH-A0HK-01A TCGA-AO-A0JF-01A # # 1: TP53 -0.9326667 - -0.4386667 -0.5056667 - -0.1083333 # # 2: BRCA1 -2.0526667 - -2.0460833 -2.0060833 - -2.0717500 # # 3: NFKB1 1.1812000 - 1.3117000 0.8047000 - 1.1610000 # # TCGA-BH-A0DP-01A TCGA-AO-A0JJ-01A TCGA-A8-A0A6-01A TCGA-AO-A0JM-01A # # 1: # # 2: -0.182500 -0.3288333 -0.535500 -0.06183333 -1.448833 - -0.9750000 -1.673083 - -1.78058333 # # 3: 1.490200 1.5671000 1.358000 0.62590000 # # TCGA-BH-A0HB-01A # # 1: -0.7593333 # # 2:-2.3987500 ## 3: 0.8151000

0.4.4映射到基因位置的关节CNA状态表

第三，读取、收集CNA数据，并使用基因组参考hg19将其映射到人类基因位置。每个基因将被评估为CNA-gain (1)， CNA-loss(-1)和neutral(0)。可以通过设置Threshold来调整CNA增益或损失的水平。

这里我们为增益事件设置阈值2.4，为损失事件设置阈值1.6。

My_cna_reader <- function(cna_filepath) {cna <- data. My_cna_reader <- function(cna_filepath) {table::fread(cna_filepath, sep = '\t'， header = TRUE) cna[， .(染色体，Start, End, Segment_Mean)]} gain_loss = log2(c(2.4, 1.6)) - 1 gene_cna <- create_gene_cna(sample_desc, gain_threshold = gain_loss[1]， loss_threshold = gain_loss[2]， read_fun = my_cna_reader, progress =FALSE) gene_cna[GENE %in% c('TP53'， 'BRCA1'， 'NFKB1')， 1:10, with=FALSE]
# #基因TCGA-AO-A0JL-01A TCGA-BH-A0HF-01A TCGA-BH-A0HK-01A TCGA-AO-A0JF-01A # # 1: BRCA1 1 1 0 0 # # 2: NFKB1 0 1 0 0 # # 3: TP53 0 0 1 0 # # TCGA-BH-A0DP-01A TCGA-AO-A0JJ-01A TCGA-A8-A0A6-01A TCGA-AO-A0JM-01A # # 1: 0 0 0 1 # # 2: 0 0 0 0 # # 3: 1 0 0 1 # # TCGA-BH-A0HB-01A # # 1: 1 # # 2: 0 # # 3: 1

0.4.4.1并行化

对于性能问题，可以启用并行化来提高流程。

#更改4以匹配一个人的总CPU核心doMC::registerDoMC(cores = 4) gene_cna <- fast_gene_cna (sample_desc, gain_loss[[1]]， gain_loss[[2]]， parallel = TRUE)

0.4.4.2直接从数据中读取基因信息

如果一个人的CNA数据已经包含基因信息，请尝试direct_gene_cna．

0.4.5cna驱动的差异表达基因鉴定

最后,运行find_cna_driven_gene从两种加工过的GE样品中鉴定由CNAs驱动的差异表达基因gene_exp和CNA数据gene_cna从前面的步骤中得到的。

给出拷贝数改变样本(增益或损失)比例的阈值，以选择cna驱动基因。

cna_driven_genes <- find_cna_driven_gene(gene_cna, gene_exp, gain_prop = 0.15, loss_prop = 0.15, progress = FALSE, parallel = FALSE) head(cna_driven_genes$gain_driven)
## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## #0.0625 2.6232500 0.6852500 -0.670500 ## 3: 0.1250 -0.7475000 -1.5413333 -1.505833 ## 4: 0.0000 1.2829792 0.4458750 NaN ## 5: 0.1250 0.2345417 1.5736818 2.524063 ## 6: 0.0625 1.5993333 0.6277917 -0.069000 ## vs_rest_exp_diff ## 1: 0.7096528 ## 2: 2.0509792 ## 3: 0.7883718 ## 4: 0.8371042 ## 5: -1.4853526 ## 6: 1.0251410
头(cna_driven_genes loss_driven美元)
##基因p_value fdr gain_sample_prop normal_sample_prop ## 1: POLK 8.417198e-10 3.447685e-06 0.1250 0.6875 ## 2: EZH1 4.227489e-07 7.088258e-04 0.0625 0.7500 ## 3: SOD2 5.191596e-07 8.2488258e -04 0.1250 0.6875 ## 5: NBR1 1.053985e-06 8.245376e-04 0.0625 0.5000 ## 6: ZNF554 1.207819e-06 8.245376e-04 0.1250 0.6875 ## loss_sample_prop gain_exp_mean normal_exp_mean loss_exp_mean ## 1: 0.1875 -0.403300 -0.1174404 -1.2971667 ## 2:0.4375 NaN 0.2341667 -0.4589762 ## 3: 0.1875 0.066250 -0.2052812 -1.0430417 ## 4: 0.1875 0.364375 -0.0490000 -1.0135833 ## 5: 0.4375 0.644500 0.8796875 -0.3331786 ## 6: 0.1875 0.950625 0.6446591 -0.2225833 ## vs_rest_diff ## 1: -1.1357479 ## 2: -0.6931429 ## 3: -0.8586474 ## 4: -1.0281795 ## 5: -1.1867341 ## 6: -0.9143141
头(cna_driven_genes两美元)
##基因gain_p_value gain_fdr loss_p_value loss_fdr gain_sample_prop ## 1: AHNAK2 0.0001146049 0.01676318 0.103289129 0.32975236 0.1875 ## 2: NTF3 0.0002403053 0.02343285 0.559826577 0.76690624 0.1875 ## 3: ABR 0.0011628790 0.05807826 0.001161531 0.04115464 0.1875 ## 4: RFX3 0.0011221713 0.05807826 0.107495688 0.33610865 0.1875 ## 5: IDH3B 0.0014698661 0.06051011 0.658967508 0.83357965 0.1875 ## 6:FNTB 0.0016109439 0.06315896 0.06396669 0.1875 ## ## normal_sample_prop loss_sample_prop gain_exp_mean normal_exp_mean ## # 1: 0.625 0.1875 0.8897222 -0.26045227 ## # 2: 0.625 0.1875 -2.8278333 -1.33816667 ## # 3: 0.375 0.4375 1.0030000 0.641377500 ## 4: 0.625 0.1875 0.7407222 -0.05018333 ## 5: 0.625 0.1875 -0.5640417 -0.98650000 ## 6: 0.625 0.1875 0.7234333 0.33618000 ## # loss_exp_mean gain_vs_rest_exp_diff ## 1: -1.08888889 1.3413522 -1.0938615 ## 2:-1.20477778 -1.5204487 0.4771581 ## 3: -0.02914286 0.7226731 -0.7910595 ## 4: -1.00261111 1.0106966 -1.1349444 ## 5: -1.03441667 0.4335160 -0.1454071 ## 6: -0.18523333 0.5075795 -0.6107795

例如，在这个数据集中，BRCA1似乎是CNA损失驱动的，它在CNA损失样本中的表达低于CNA中性样本。

cna_driven_genes$loss_driven[基因%in% c('BRCA1')]
## GENE p_value fdr gain_sample_prop normal_sample_prop ## 1: BRCA1 0.3259034 0.5930255 0 0.5625 ## loss_sample_prop gain_exp_mean normal_exp_mean loss_exp_mean ## 1: 0.4375 NaN -1.547074 -1.773119 ## vs_rest_exp_diff ## 1: -0.226045

0．5问与答

0.5.0.1问:为什么要求使用捆绑的hg19人类基因组参考?

在发展的早期，iGC需要一个特殊的数据结构作为基因组参考，因此捆绑了一个。现在不需要这样的特殊结构，所以我们计划在即将发布的版本中放松这样的限制，用户将能够传递在Bioconductor上可用的其他参考资料。

目前修改的hg19DBNM包含UCSC Genome Browser中hg19的RefSeq转录本。选择具有NM标记ID和蛋白编码的转录本。