转录组生信分析教程
桓峰基因公众号推出转录组分析教程，有需要生信的老师可以联系我们！转录分析教程整理如下：
RNA 2. SCI文章中基于GEO的差异表达基因之 limma
RNA 3. SCI 文章中基于T CGA 差异表达基因之 DESeq2
RNA 4. SCI 文章中基于TCGA 差异表达之 edgeR
RNA 6. 差异基因表达之-- 火山图 (volcano)
RNA 7. SCI 文章中的基因表达——主成分分析 (PCA)
RNA 8. SCI文章中差异基因表达--热图 (heatmap)
RNA 12. SCI 文章中肿瘤免疫浸润计算方法之 CIBERSORT
RNA 14. SCI 文章中差异表达基因之 蛋白互作网络 (PPI)
RNA 15. SCI 文章中的融合基因之 FusionGDB2
RNA 17. SCI 文章中的筛选 Hub 基因 (Hub genes)
RNA 19. SCI 文章中无监督聚类法 （ConsensusClusterPlus）
RNA 20. SCI 文章中单样本免疫浸润分析 （ssGSEA）
RNA 22. SCI 文章中基于表达估计恶性肿瘤组织的基质细胞和免疫细胞(ESTIMATE)
RNA 23. SCI文章中表达基因模型的风险因子关联图(ggrisk)
RNA 24. SCI文章中基于TCGA的免疫浸润细胞分析 (TIMER)
RNA 25. SCI文章中估计组织浸润免疫细胞和基质细胞群的群体丰度(MCP-counter)
RNA 26. SCI文章中基于转录组数据的基因调控网络推断 (GENIE3)
RNA 27 SCI文章中转录因子结合motif富集到调控网络 (RcisTarget)
RNA 28 SCI 文章中基于RNA-seq数据反褶积揭示肿瘤免疫结构的分子和药理学 (quanTIseq)
RNA 29. SCI文章中基于TCGA的免疫浸润细胞分析 (TIMER2.0)
RNA 30. SCI文章中基于TCGA和GTEx数据挖掘神器(GEPIA2)
RNA 31. SCI文章临床蛋白质组肿瘤在线数据挖掘神器（CPTAC）
RNA 32. SCI文章临床多组学肿瘤在线数据挖掘神器（UALCAN）
RNA 33. SCI文章肿瘤在线数据挖掘神器（cBioportal）
RNA 34. SCI 文章基因蛋白表达时间趋势分析及划分聚类群 (Mfuzz)
简 介
转录组测序(RNA-seq)被广泛用于检测急性淋巴细胞白血病(ALL)的基因重排和定量基因表达，但其在鉴定CNVs方面的实用性和准确性尚未得到很好的描述。由于非整倍体亚型在ALL中的高发，通过RNA-seq推断出的CNV信息对于指导ALL的疾病分类和风险分层具有很高的信息性。在这里，我们描述了RNAseqCNV，一种从RNA-seq数据中检测大规模拷贝数变化的方法。我们使用基于归一化基因表达和小等位基因频率的模型对ALL中臂水平的CNVs进行了分类，准确率很高(总体为99.1%，非二倍体染色体臂为98.3%)，该模型在急性髓性白血病中进一步验证了出色的表现(总体为99.8%，非二倍体染色体臂为99.4%)。RNAseqCNV在调用ALL数据集中的cnv方面优于其他基于RNA-seq的算法，特别是在具有高比例cnv的样本中。CNV与基于DNA的CNV结果高度一致，比传统的基于细胞遗传学的核型更可靠。RNAseqCNV提供了一种在没有DNA分析的情况下稳健识别拷贝数变化的方法，进一步增强了RNA-seq对ALL亚型的分类能力。
软件包安装
# install devtools
install.packages("devtools")
# install RNAseqCNV package
devtools::install_github(repo = "honzee/RNAseqCNV")数据读取
软件包和Shiny应用程序都接收到相同的输入。每个基因读取计数和SNV突变等位基因频率(MAF)被用来产生结果。因此，每个样本需要两种类型的信息。可以在包目录中找到输入文件的示例。使用如下命令定位软件包目录:
# Examples of basic function calls:
library(RNAseqCNV)
system.file("extdata", package = "RNAseqCNV")
## [1] "D:/R-4.2.1/library/RNAseqCNV/extdata"
list.files(system.file("extdata", package = "RNAseqCNV"))
## [1] "config_mock.txt" "diploid.HTSeq"
## [3] "diploid.snv" "high_hyperdiploid.HTSeq"
## [5] "high_hyperdiploid.snv" "low_hypodiploid.HTSeq"
## [7] "low_hypodiploid.snv" "metadata_mock"
## [9] "near_haploid.HTSeq" "near_haploid.snv"要运行软件包或Shiny应用程序，需要正确格式的配置文件和metadata文件的路径:
config 文件
配置文件定义了输出目录(out_dir)、读计数文件目录(count_dir)和SNV文件目录(snv_dir)。更改相应的文件目录，但保持关键字相同，如下例所示:
out_dir = "/Path/to/output_dir"
count_dir = "/Path/to/dir/with/count_files"
snv_dir = "/Path/to/dir/with/vcf_files"比如我们这里使用的是软件包里面的例子：
将这个文件copy到自己的分析目录，否则会保存，因为输出目录就是当前目录。
out_dir = "./"
count_dir = system.file("extdata", package = "RNAseqCNV")
snv_dir = system.file("extdata", package = "RNAseqCNV")metadata 文件
metadata文件包含三列(以逗号/tab/空格分隔):
1.样品标识;
2.计算文件;
3.Vcf /custom文件。
注：不接受标题行，列的顺序必须与下面的示例相同:
diploid diploid.HTSeq diploid.snv
high_hyperdiploid high_hyperdiploid.HTSeq high_hyperdiploid.snv
low_hypodiploid low_hypodiploid.HTSeq low_hypodiploid.snv
near_haploid near_haploid.HTSeq near_haploid.snv再来详细说明以下这三种文件的格式：
Read count file
HTSeq-count或类似读计数软件的输出。表必须有两列:
1.集合基因id;
2.read counts;
注：表不应该有标题。
ENSG00000000005 0
ENSG00000000419 94
ENSG00000000457 128
ENSG00000000460 133
ENSG00000000938 171
ENSG00000000971 5SNV information
VCF文件可以通过运行GATK流程获得的vcf文件：
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample_name
1 14792 . G A 156.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=-1.442;ClippingRankSum=0.000;DP=27;ExcessHet=3.0103;FS=1.690;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=5.81;ReadPosRankSum=-0.169;SOR=1.124 GT:AD:DP:GQ:PL 0/1:19,8:27:99:185,0,644
1 14907 . A G 126.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=-0.015;ClippingRankSum=0.000;DP=10;ExcessHet=3.0103;FS=2.808;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=14.09;ReadPosRankSum=0.751;SOR=1.170 GT:AD:DP:GQ:PL 0/1:3,6:9:64:155,0,64
1 14930 . A G 161.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.751;ClippingRankSum=0.000;DP=10;ExcessHet=3.0103;FS=2.808;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=17.97;ReadPosRankSum=2.515;SOR=1.170 GT:AD:DP:GQ:PL 0/1:3,6:9:64:190,0,64Custom tabular data
该表有四个必需的列:
chr: 染色体;
start: SNV的位置;
depth: SNV位点的深度;
maf:突变等位基因频率，可以计算为突变等位基因的深度(与参考基因组相比)除以该位点的总体 reads的深度。
head名称应遵循以下格式:
chr start depth maf
1 14599 40 0.5
1 14604 9 0.3333
1 14610 10 0.25实操
基因组选择
可设置参数：genome_version = "hg38" or "hg19" 。
# example of wrapper function (DO NOT RUN)
RNAseqCNV_wrapper(config = "config_mock.txt", metadata = "metadata_mock", snv_format = "custom",
genome_version = "hg38")Arm-level figures
Arm-level figures的绘制显著增加了每个样本的运行时间。您可以通过arm_lvl = FALSE减少运行时间:
RNAseqCNV_wrapper(config = "config_mock.txt", metadata = "metadata_mock", snv_format = "custom",
arm_lvl = FALSE)Estimation labels
可以使用以下方法从图中删除CNV估计标签:
RNAseqCNV_wrapper(config = "config_mock.txt", metadata = "metadata_mock", snv_format = "custom",
estimate_lab = FALSE)Diploid adjustment
有些样本可能具有较高比例的CNVs染色体，如下图所示:
考虑到RNA-seq数据的相对性质，在这些样本中，稳定表达的基因在二倍体染色体上的归一化表达不会以零为中心。为了解决这个问题，软件包包括一个随机森林模型，它将染色体分类为二倍体或非二倍体。根据这些信息，图形的居中方式与估计为二倍体的染色体居中于零的方式相同:
此功能在默认情况下是打开的，而不需要矫正的话就adjust = FALSE:
RNAseqCNV_wrapper(config = "config_mock.txt", metadata = "metadata_mock", snv_format = "custom",
adjust = FALSE)Analysis mode
样品可以分析两种模式:批量分析或每个样品分析。每个样本分析内置标准样本的基因表达规范化和定心是默认的。在批量分析中，输入样本将用于归一化和基因表达中心(log2倍变化计算)。在这种模式下，至少需要提供20个没有大量大规模CNVs的样本。为获得最佳结果，样品应具有相同的癌症类型和文库制备。要使用这种分析模式，请在RNAseqCNV_wrapper函数中使用batch = TRUE。
RNAseqCNV_wrapper(config = "config_mock.txt", metadata = "metadata_mock", snv_format = "custom",
batch = TRUE)在每个样品分析中，每个样品将根据内置(或用户提供)标准样品进行分析。内建标准包含来自40个无大规模CNVs的ALL样本的基因表达数据。使用内置标准示例是默认方法。也可以提供标准样品作为输入。标准样本必须包含在metadata表中，并且它们的样本名称必须在RNAseqCNV_wrapper函数中以字符向量格式(standard_samples参数)提供。
RNAseqCNV_wrapper(config = "config_mock.txt", metadata = "metadata_mock", snv_format = "custom",
batch = FALSE, standard_samples = c("standard_sample_1,", "standard_sample_2,",
"standard_sample_3,"))Gene weights generation
除CNVs外，影响基因表达水平的因素还有很多。因此，RNAseqCNV基于精心策划的CNV数据集为每个基因分配一个权重值，以利用每个基因在预测CNV中的重要性。默认的内置基因权重基于426个B-ALL样本，源于基因在表达水平上的cnv相关性和方差。要从输入数据生成权重，参数generate_weights应设置为TRUE。然而，这些权重只考虑了样本间基因表达的差异。
RNAseqCNV_wrapper(config = "config_mock.txt", metadata = "metadata_mock", snv_format = "custom",
generate_weights = TRUE)也可以根据以下格式提供自定义权重矩阵:
ENSG chromosome_name weight
ENSG00000000419 20 121.923504
ENSG00000001167 6 8.474673
ENSG00000001497 X 442.760032R函数get_weights():
file.path(find.package("RNAseqCNV"), "R", "user_defined_analysis", "get_weights")CNV matrix
可选地，估计的CNVs也可以以矩阵格式输出，其中样本按行排序，染色体臂按列排序。
RNAseqCNV_wrapper(config = "config_mock.txt", metadata = "metadata_mock", snv_format = "custom",
CNV_matrix = TRUE)得到的文件格式如下：
sample 1p 1q 2q 2p 3p 3q 4p 4q 5q 5p 6p 6q
sample_1 0 0 0 0 0 0 0 0 0 0 0 0
sample_2 0 0 0 0 0 0 0 0 0 0 0 0
sample_3 0 0 0 0 0 0 1 1 0 0 1 1
sample_4 0 0 0 0 0 0 0 0 0 0 1 1另一方面，Shiny应用程序可以更轻松地浏览和检查结果。
launchApp()结果说明
软件包和Shiny应用程序的输出(图形和表格)都将保存在配置文件中指定的输出目录中。每个样本会有一个输出目录，里面都是图片，另外还有三张表格：
Main figure
主体图形由两个面板组成。上图显示了每条染色体的基因表达水平。Y轴表示基因表达相对于参比样品的log2倍变化;x轴显示1-22-X的23条染色体(Y染色体除外)。X轴也表示基因在每条染色体上的位置。对于每条染色体，根据归一化基因表达的分布绘制加权箱线图(1/4、1/2和3/4分位数)。对于每条染色体，使用随机福雷斯特算法估计CNVs，并在每条染色体上标记结果。低置信度(质量)的CNV用红色突出显示。下图显示了每条染色体的MAF密度图。请注意，仅使用杂合snv的MAF (MAF在0.05 ~ 0.9之间)来测定cnv。峰距，测量x轴上MAF密度图中两个最高峰之间的距离，也标记在顶部。
Arm-level figures
用户可以选择生成Arm-level的CNV数据：
RNAseqCNV_wrapper(config = "config_mock.txt", metadata = "metadata_mock", snv_format = "custom", arm_lvl = TRUE)
Estimation table
估计的性别、臂位变化和染色体数目保存在输出目录的两个表中。estimation_table.tsv存储RNAseqCNV模型的输出。manual_an_table.tsv通过Shiny应用程序存储用户手动策划的结果。
Sample gender chrom_n alterations
SJHYPER141_D male 59 1q+, 4+, 5+, 6+, 7+, ?8+, ?8+, 10+, 12+, 14+, 14+, ?16q, 17+, ?18+, ?18+, 21+, 21+, 22+, X+
SJALL015971_D1 female 28 1-, 2-, 3-, 4-, 5-, 7-, 8-, 9-, 11-, 12-, 13-, 14-, 15-, 16-, 17-, 19-, 20-, 22-
SJALL049672_D1 male 34 2-, 3-, 4-, 5-, 7-, 9-, 13-, 15-, 16-, 17-, 20-, ?21, 22-
SJALL015927_D1 female 52 ?6p+, 6q+, 10+, 14+, 17q+, 18+, 21+, 21+, X+Reference
Bařinka J, Hu Z, Wang L, et al. RNAseqCNV: analysis of large-scale copy number variations from RNA-seq data. Leukemia. 2022;36(6):1492-1498. doi:10.1038/s41375-022-01547-8
这个工具是不是非常好用？
学会了轻松拿捏数据分析，顿时觉得自己也行了，关注桓峰基因公众号，轻松学生信，高效写文章！
桓峰基因，铸造成功的您！
未来桓峰基因公众号将不间断的推出单细胞系列生信分析教程，
敬请期待！！
桓峰基因官网正式上线，请大家多多关注，还有很多不足之处，大家多多指正！
http://www.kyohogene.com/
桓峰基因和投必得合作，文章润色优惠85折，需要文章润色的老师可以直接到网站输入领取桓峰基因专属优惠券码：KYOHOGENE，然后上传，付款时选择桓峰基因优惠券即可享受85折优惠哦！https://www.topeditsci.com/
扫一扫
2777
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
