物联网app开发 单个基因的拷贝数变异频率（gain or loss CNV frequency）棒棒糖图
发布日期：2024-11-06 09:03 点击次数：181
💡专注R话语在🩺生物医学中的使用
设为“星标”，精彩可以过
众人在一些触及拷贝数变异（copy number variation，cnv）的SCI著述中可能会遭遇以下图片：
DOI:10.1007/s00011-023-01772-6
图片
https://doi.org/10.1186/s12943-020-01170-0图片
这个图形展示了单个基因的拷贝数变化。如若你也曾有了数据，绘图照旧很直率的，等于一个直率的点和线的组结伴料，有多种容貌可以竣事，比如我之前就成心先容过：
R话语棒棒糖图R话语画瑕玷线的5种期间今上帝要记载下若何通过拷贝数变异数据，取得绘图的数据。咱们就以从TCGA下载的TCGA-COAD的拷贝数变异数据为例。
加载数据径直从TCGA下载的拷贝数变异数据是不行的，需要使用GISTIC2.0这个软件对CNV数据进行分析，取得以下成果才行：
图片
GISTIC2.0的输出这个经由对于入门者照旧很复杂的，光是装配gistic这个软件测度就会卡住许多东说念主。可是这个经由网上有极端多的教程，咱们就未几作念先容了。这里给众人提供一个装配GISTIC软件的快速期间，可以1行代码装配，详见：王诗翔大佬的github：Install GISTIC2 by one line code
除此以外，这个分析还可以通过网页器用完成，一个是gene pattern，另一个是hiplot（不知说念当今还能不可用这个功能）。
GISTIC2.0的输出成果中有多个文献都可以用来绘制上头的图，可是它们是有区别的，其中部分红果的证据参考简书（网址：https://www.jianshu.com/p/5de9922e6b8b）：
all_data_by_genes.txt：基因在不相通本中具体的拷贝数数值all_thresholded.by_genes.txt：基因在不相通本中拷贝数数值冲破化后的成果，-2代表缺失两个拷贝，-1代表缺失一个拷贝，0代表拷贝数正常，1代表加多一个拷贝，2代表扩增两个拷贝focal_data_by_genes.txt：基因在不相通本中具体的拷贝数数值（只辩论 focal events）broad_data_by_genes.txt：基因在不相通本中具体的拷贝数数值（只辩论 arm events）all_lesions.conf_90.txt：识别到的拷贝数扩增和缺失的Peak区域amp_genes.conf_90.txt：识别到的拷贝数扩增的Peak区域及区域内触及到的基因del_genes.conf_90.txt：识别到的拷贝数缺失的Peak区域及区域内触及到的基因broad_significance_results.txt：权贵发生拷贝数变异的broad区域broad_values_by_arm.txt：染色体臂在样本中的拷贝数数值scores.gistic：该算法的打分红果，可导入IGV进行可视化其中前4个文献都可以用来绘图，可是凭证它们的具体含义我最终遴选了all_thresholded.by_genes.txt这个文献。
最初咱们详情是遴选扫数的样本，不可只辩论focal或者arm，对于这两种东西的证据，参考一个博客的证据：
在通盘基因组上，最常见的体细胞拷贝数变异（SCNA）是短的染色体变异（focal），和果然和染色体臂或整条染色体等长的染色体变异（arm-level）。arm-level的染色体变异发生的占比概况是focal的30倍，且果然扫数的癌症类型中都是这么。合计这两者发生概率不同，意味着是分手发生的。https://taozyblog.com.cn/post/2022-3-16_mechanism_of_cnv/
其次是使用具体的拷贝数数值照旧冲破化之后的数据，我都读取进来对比了一下：
1. 浦项铁人俱乐部成立于1973年，球队历史曾获得5次韩K联赛冠军，4次韩国杯冠军，2次韩国联赛杯冠军，物联网app开发1次亚冠联赛冠军，以及在96/97/和97/98连续获得亚冠前身亚洲俱乐部锦标赛冠军等诸多赛事荣誉。
# 具体的拷贝数数值df1 <- data.table::fread("G:/tcga/TCGA_COAD_results/all_data_by_genes.txt", data.table = F)dim(df1)## [1] 25988 492
df1[1:4,1:4]
## Gene Symbol Gene ID Cytoband TCGA-DM-A285-01A-11D-A16U-01## 1 ACAP3 116983 1p36.33 -0.263## 2 ACTRT2 140625 1p36.32 -0.263## 3 AGRN 375790 1p36.33 -0.263## 4 ANKRD65 441869 1p36.33 -0.263
# 冲破化后的数据df2 <- data.table::fread("G:/tcga/TCGA_COAD_results/all_thresholded.by_genes.txt", data.table = F)dim(df2)## [1] 25988 492
df2[1:4,1:4]
## Gene Symbol Locus ID Cytoband TCGA-DM-A285-01A-11D-A16U-01## 1 ACAP3 116983 1p36.33 -1## 2 ACTRT2 140625 1p36.32 -1## 3 AGRN 375790 1p36.33 -1## 4 ANKRD65 441869 1p36.33 -1
可以看到它们的数据结构是裕如一样的，仅仅暗示期间不一样。如若你遴选具体的拷贝数数值，那么你需要我方设定一个阈值（0.3，0.4，0.5等），到底这个数值大于若干是gain（也等于duplication），小于若干是loss（也等于deletion），可是如若使用冲破化之后的数据，那么凭证它的意念念，大于0的等于gain，小于0的等于loss，不需要我方设定阈值了，为了通俗，我遴选这个冲破化之后的数据。
数据整理这个df2是长这么的：
图片
第一列是基因名字，第2和第3列咱们绘图用不到，可以径直删掉，背面的列都是样真名，其中是用数字暗示拷贝数变异，是以咱们对它进行极幼年小的整理：
去掉第2和第3列把-1和-2酿成loss，把1和2酿成gain，0酿成neutraldf2 <- df2[,-c(2,3)]rownames(df2) <- df2[,1]df2 <- df2[,-1]df2[df2>0] <- "gain"df2[df2<0] <- "loss"df2[df2 == 0] <- "neutral"
然后某个基因的CNV frequency就很好运筹帷幄了，比如第一个基因它的gain的频率，等于出现gain的样本数 / 样本总和，它的loss的频率等于出现loss的样本数 / 样本总和，先直率看下出现gain/loss/neutral的样本数目，以第一个基因为例：
table(t(df2)[,1])
## ## gain loss neutral ## 18 162 309
凭证这个念念路咱们就可以很平常的运筹帷幄出扫数基因的gain/loss/neutral的样本数：
df3 <- as.data.frame(t(df2))df3 <- apply(df3, 2, table)df3 <- as.data.frame(do.call(rbind,df3))
## Warning in (function (..., deparse.level = 1) : number of columns of result is## not a multiple of vector length (arg 24509)
df3$gene <- rownames(df3)head(df3)
## gain loss neutral gene## ACAP3 18 162 309 ACAP3## ACTRT2 18 162 309 ACTRT2## AGRN 18 162 309 AGRN## ANKRD65 18 162 309 ANKRD65## ATAD3A 18 162 309 ATAD3A## ATAD3B 18 162 309 ATAD3B
这个等于咱们的绘图数据了（我这里并莫得除以样本总和，众人有需要的可以我方运筹帷幄），有了这个数据绘图等于径直ggplot2就好了。
绘图咱们决然取30个基因，你可以凭证我方的情况遴选我方感意思的基因，比如免疫联系的基因、铜牺牲联系的基因等。
set.seed(1567)markers <- sample(df3$gene, 30)markers
## [1] "KRT10" "LAMA1" "TADA3" "MIR663A|chr20"## [5] "LOC100129138" "hsa-mir-886" "DSC2" "KRTAP22-2" ## [9] "MICU3" "LOC390705" "ZYG11B" "LRRC75A" ## [13] "FAM86B1" "TRIM52-AS1" "LOC102723604" "MIR4304" ## [17] "RNF139-AS1" "LINC00620" "ISCA1" "FBP1" ## [21] "TRIAP1" "APOA1-AS" "MLN" "ST6GALNAC2" ## [25] "OR4K17" "MIR518C" "GABPB2" "ZFP28" ## [29] "MEIOC" "CSRP1"
df4 <- df3[df3$gene %in% markers,]df4
## gain loss neutral gene## ZYG11B 12 130 347 ZYG11B## LOC100129138 12 151 326 LOC100129138## GABPB2 113 37 339 GABPB2## CSRP1 115 42 332 CSRP1## TADA3 39 75 375 TADA3## LINC00620 39 76 374 LINC00620## hsa-mir-886 49 99 341 hsa-mir-886## TRIM52-AS1 57 94 338 TRIM52-AS1## MLN 104 42 343 MLN## FAM86B1 92 217 180 FAM86B1## MICU3 96 216 177 MICU3## RNF139-AS1 283 12 194 RNF139-AS1## ISCA1 73 62 354 ISCA1## FBP1 76 64 349 FBP1## APOA1-AS 61 73 355 APOA1-AS## TRIAP1 101 55 333 TRIAP1## MIR4304 102 55 332 MIR4304## OR4K17 42 149 298 OR4K17## LOC102723604 41 149 299 LOC102723604## LOC390705 131 16 342 LOC390705## LRRC75A 13 276 200 LRRC75A## KRT10 120 52 317 KRT10## MEIOC 108 59 322 MEIOC## ST6GALNAC2 112 62 315 ST6GALNAC2## LAMA1 25 279 185 LAMA1## DSC2 28 288 173 DSC2## MIR518C 97 56 336 MIR518C## ZFP28 102 53 334 ZFP28## MIR663A|chr20 290 29 170 MIR663A|chr20## KRTAP22-2 30 141 318 KRTAP22-2
绘图是最直率的一步：
library(ggplot2)ggplot(df4)+ geom_linerange(aes(x=gene,ymin=0,ymax=ifelse(gain>loss,gain,loss)), linewidth=2,color="grey")+ geom_point(aes(x=gene, y=gain), color="red",size=4)+ geom_point(aes(x=gene, y=loss), color="blue",size=4)+ labs(x=NULL,y="CNV frequency")+ theme_bw()+ theme(axis.text.x = element_text(angle = 45,hjust = 1))
图片
你还可以先排个序再画，这么看起来更整王人极少物联网app开发，我这里就不演示了。
本站仅提供存储奇迹，扫数本色均由用户发布，如发现存害或侵权本色，请点击举报。