RNA-Seq 数据集、比对和标准化

RNA-Seq 数据集、比对和标准化|玉米中的元基因调控网络突出了功能上相关的调控相互作用。

RNA-Seq 表达分析代码和数据

该仓库是一个公开可用 RNA-Seq 数据集的集合（主要是玉米数据），提供了系统分析这些数据的代码/流程，以及质量控制（QC）和总结性数据输出。重点是大规模的 Illumina RNA-Seq 实验（包括多种组织/发育阶段、多个自交系/杂交系），但也包括使用其他测序技术（如 3’ RNA-Seq 等）进行的实验。

原始的测序读数从 NCBI Sequence Read Archive (SRA) 下载，使用 Trim Galore 或 fastp 进行修剪，随后使用 Hisat2 或 STAR 将其比对到玉米 B73 AGP_v4 基因组。唯一比对的读取被分配给 46,117 个参考基因模型（Ensembl Plants v37），并使用 featureCounts 进行计数。原始读取计数随后使用 TMM 标准化方法进行标准化，得到 CPM（每百万读取的计数值），然后通过基因 CDS 长度进一步标准化，得到 FPKM（每千碱基的外显子每百万读取的片段数）值。层次聚类和主成分分析（PCA）被用于探索样本聚类模式。

收集的数据集列表：

yid	author	year	source	accession	study	genotype	tissue	n	ASE	stress	RIL	Run
ca20a3	Anderson	2020	local	sp068a				59	T	T		C
cp12a2	Bolduc	2012	sra	PRJNA168086								C
cp12b2	Morohashi	2012	sra	PRJNA167802 PRJNA167803								C
cp14g2	Eveland	2014	sra									C
cp15a2	Pautler	2015	sra									C
cp15b2	Li	2015	sra									C
cp16a2	Yang	2016	sra									C
cp16b2	Gontarek	2016	sra	PRJNA260183								C
cp18a2	Zhan	2018	sra									C
cp18b2	Li	2018	sra									C
cp19c2	Dong	2019	sra		Zm00001d033673 Zm00001d028129							C
rn10a	Li	2010	sra	PRJNA79627	leaf 6 zones	B73	leaf	6				C
rn11a	Davidson	2011	sra	PRJNA80041	reproductive 12	B73	reproductive	12				C
rn13a	Li	2013	sra	PRJNA179160	eQTL	105 RILs + 2 parents	SAM	107		T		C
rn13b	Liu & Yu	2013	sra	PRJNA179196 PRJNA253977	leaf time-series 13t + 9t	B73	early and late leaf	22				C
rn13c	Eichten	2013	sra	PRJNA173886	diverse inbreds	62 inbreds	seedling_leaf3	62				C
rn14f	Makarevitch	2014	sra	PRJNA244661	stress	B73, B37, Oh43	seedling	27	T	T		C
rn15d	Baute	2015	sra	PRJEB9918	B73 x H99 RIL	104 RILs + 2 parents	leaf_4	106		T		C
rn16b	Stelpflug	2016	sra	PRJNA171684 SRP010680	atlas	B73	dev atlas	94				C
rn16c	Walley	2016	sra	PRJNA217053	proteome	B73	23 tissues	23				C
rn17b	Waters	2017	local	sp033	cold and heat stress	B, M, P, BxM, BxP	seedling_leaf3	15	T	T		C
rn17c	Marcon	2017	sra	PRJNA284670	drought stress	B, M, BxM	root	8	T	T		C
rn17d	Zhang	2017	sra	PRJNA343268 PRJNA344653	cold stress	maize, sorghum	10d whole plant	6		T		C
rn17e	Anderson	2017	local	sp060	4 stresses	B, M, P, W, O and 4 others	leaf, root	50		T		C
rn18g	Zhou	2018	local	sp052	B73 Mo17 atlas	B, M, BxM	23 tissues	90	T			C
rn18i	Patrick	2020	local	ch001	gene fusion	B, P, W	10 tissues	30				C
rn18j	Cao	2018	sra	PRJNA477643	drought stress	Yu882	leaf	6	T	T		C
rn18k	Meng	2018	sra	PRJNA377057	sequential imprinting	BxM, MxB	embryo					C
rn19a	Mazaheri	2019	sra	PRJNA437324	widiv biomass	453 inbreds	whole seedlings	453				C
rn19d	He	2019	sra	PRJNA520822	heat stress	B73	7 tisues	21	T	T		C
rn19f	Anderson	2019	local	rn19f	submergence	B73	leaf, coleoptile, root	18	T	T		C
rn19g	Hendron	2019	sra	PRJEB32558	light exposure		leaf_2	10				C
rn19h	Chang	2019	sra	PRJNA450300	comparative leaf GRN	B73	leaf	12				C
rn20a	Zhou	2020	local	sp069	cold and heat stress time series	30 genotypes	leaf	291	T	T		C
rn20a2	Zhou	2020	local	rn20a2	BRB-Seq				T			C
rn20b	Zhou	2020	local	sp065	heterosis	4 inbred + 6 hybrids	3 tissues		T			C
rn20b2		local	sp068a	RIL	B, M, BxM, MxB, 4 RILs	leaf		T	T		C
rn20c		2020	sra	PRJEB36014	NAM25	25 NAM parents	10 tissues	250				C
rn20d	Li	2019	sra		biomAP	202 inbreds + 96 hybrids	6 tissues	636				C
rn20d3	Li	2019	local

输出文件的详细说明

样本列表 / 元数据表：`01.meta.tsv`

SampleID：样本ID
Tissue：组织类型
Genotype：基因型
Treatment：处理方法
Replicate：重复编号
paired：是否为双端（paired-end）测序
spots：读取数（单端）或对数（双端）
avgLength：平均读长

结果文件，包括修剪和映射QC统计、原始读数计数及标准化的CPM / FPKM表：

01.rds

可通过 x = readRDS("01.rds") 在R中加载，包含以下数据框：
- th：样本列表 / 元数据表，与 01.meta.tsv 相同
- trimming：修剪统计
  - sid：样本ID
  - passed_filter_reads：通过过滤的读取数
  - low_quality_reads：低质量的读取数
  - too_many_N_reads：包含过多N的读取数
  - too_short_reads：过短的读取数
  - too_long_reads：过长的读取数
- bamstat：映射统计
  - sid：样本ID
  - pair：双端对
    - pair_bad：失败QC的对数
    - pair_dup：重复的对数
    - pair_map：映射的对数（两个端）
    - pair_orphan：仅一个端被映射的对
    - pair_unmap：未映射的对
  - unpair：单端（单端读取或一个端未通过QC的对）
    - unpair_bad：失败QC的单端
    - unpair_dup：重复的单端
    - unpair_map：映射的单端
    - unpair_unmap：未映射的单端
  - pair_map_hq，pair_orphan_hq，unpair_map_hq：高质量（即唯一）映射的对数/读取数
  - pair_map0，pair_orphan0，unpair_map0：映射并且没有错配的对数/读取数
  - pair_map_hq0，pair_orphan_hq0，unpair_map_hq0：高质量（即唯一）且没有错配的对数/读取数
- fcnt：来自 featureCounts 的原始读取计数
  - gid：基因ID（AGP_v4，Ensembl Plants v37，总共46,117个）
  - SampleID：样本ID
  - ReadCount：原始读取计数
- salmon：来自salmon的原始读取计数和标准化的TPM（每百万转录本数）
  - gid：基因ID（AGP_v4，Ensembl Plants v37，总共46,117个）
  - SampleID：样本ID
  - ReadCount：原始读取计数
  - TPM：salmon标准化的转录本每百万数值（TPM）
- salmon_tx：来自salmon的原始读取计数和标准化的TPM（转录本级别）
  - tid：转录本ID（AGP_v4，Ensembl Plants v37）
  - SampleID：样本ID
  - ReadCount：原始读取计数
  - TPM：salmon标准化的转录本每百万数值（TPM）
- tl：文库统计
  - SampleID：样本ID
  - libSize：文库大小
  - sizeFactor：DESeq2文库大小因子
  - normFactor：edgeR文库标准化因子
- tm：标准化表达表
  - gid：基因ID（AGP_v4，Ensembl Plants v37，总共46,117个）
  - SampleID：样本ID
  - ReadCount：原始读取计数
  - nRC：标准化读取计数（nRC = ReadCount / sizeFactor）
  - rCPM：原始CPM（每个样本/文库的总和为1,000,000）
  - rFPKM：原始FPKM（使用rCPM和基因外显子长度计算）
  - rTPM：原始TPM（每个样本/文库的总和为1,000,000）
  - CPM：通过edgeR计算的CPM（CPM = rCPM / normFactor）
  - FPKM：通过CPM和基因外显子长度计算的FPKM
  - TPM：标准化的TPM（TPM = rTPM / normFactor）
- th_m：重复合并的样本列表 / 元数据表
- tm_m：重复合并的表达表
- ase_gene：基因级别的等位基因特异性读取计数
  - sid：样本ID
  - gid：基因ID
  - allele1，allele2：每个等位基因的特异性读取计数。例如，在Mo17xB73的情况下，allele1表示Mo17（第一）等位基因的计数，而allele2表示Mo17（第二）等位基因的计数
- ase_snp：SNP级别的等位基因特异性读取计数
  - sid：样本ID
  - chr，pos，ref，alt：SNP信息
  - gt：样本在此位点的基因型（0|1或1|0）
    - 在1|0的情况下，等位基因1（母源等位基因）在alt状态，而等位基因2（父源等位基因）在ref状态
  - allele1，allele2：母源（第一）等位基因和父源（第二）等位基因的读取计数

biomap

biomAP mRNA-Seq 数据集分析

方法

结果

原始样本元数据表
修正后的样本元数据表：
- 样本 bm252 从 Root 修正为 Leaf
- MSI 路径：/home/springer/zhoux379/projects/rnaseq/data/05_read_list/me99c.c.tsv
中间文件都存放在 MSI 临时空间目录下：
/scratch.global/zhoux379/rnaseq/me99c/，包括以下子目录：
- 10_fastq, 15_trim：原始和修剪后的 fastq 文件
- 21_star, 22_bam：原始 BAM 文件和坐标排序后的 BAM 文件
- 31_featurecounts：读取计数表
- 31_mmquant：原始读取计数
- 33_ase：等位基因计数分析的中间文件

QC (质量控制)：

[读取修剪和映射统计]
PCA 图
t-SNE 图：
- 相较于 PCA 图，t-SNE 图显示了更好的组织分离效果
- 叶片样本被分为两个不同的簇
映射统计表：
包含每个样本的修剪、映射和计数统计，表格列出：
- MSI 路径：/home/springer/zhoux379/projects/rnaseq/data/raw_output/bamstats.tsv
- SampleID, Tissue, Genotype, Treatment, Replicate：样本元数据
- 修剪统计：total, surviving, surviving_f, surviving_r, dropped
- 映射统计：
  - pair：读对
    - pair_bad, pair_dup：映射失败的读对或重复读对
    - pair_map：映射的读对（两端）
    - pair_orphan：只有一端映射的读对
    - pair_unmap：未映射的读对
  - unpair：单端读取（或一端失败的配对）
    - unpair_bad, unpair_dup：映射失败的单端或重复单端
    - unpair_map：映射的单端读取
    - unpair_unmap：未映射的单端读取
  - pair_map_hq, pair_orphan_hq, unpair_map_hq：高质量（即唯一）映射的读对/单端
  - pair_map0, pair_orphan0, unpair_map0：无错配的映射读对/单端
  - pair_map_hq0, pair_orphan_hq0, unpair_map_hq0：高质量并且无错配的映射读对/单端
- 读取计数统计：
  - Assigned：分配给外显子区域并计数的读取
  - Unassigned_MultiMapping, Unassigned_NoFeatures, Unassigned_Ambiguity, Unassigned_Unmapped：由于各种原因未计数的读取

ASE（等位基因特异性表达）分析：

每个样本的冲突读取比例：

大部分样本中基因的冲突读取比例非常低（通常小于 2%），异常高的冲突读取比例可能表明样本基因型标注错误。
每个样本的父本等位基因比例：

数字表示样本数（即基因数）。大多数自交系样本的父本等位基因比例为 0，而杂交样本的父本等位基因比例为 0.5，例外情况可能表明样本基因型标注错误。

R 数据文件

MSI 路径：/home/springer/zhoux379/projects/biomap/data/41_qc/10.rc.ase.rda
包含原始读取计数表、标准化表达值和等位基因特异性读取计数：
th - 样本信息表（tibble），包含以下列：
- SampleID：样本编号（如：bm001 - bm467）
- Tissue：样本组织（如：Leaf, Internode, Root 等）
- Genotype：基因型（如：B73, Mo17xPH207 等）
- Treatment：处理方式（如：replicate 1 或 2）
- inbred：是否为自交亲本（TRUE 或 FALSE）
- sizeFactor, libSize：使用 DESeq2 中位数对数比率法计算的库大小和标准化因子，用于库大小的校正
- normFactor：使用 edgeR 中的 TMM 方法计算的标准化因子，不进行库大小校正
tm - biomap 表达数据表（tibble）
- gid：基因 ID（AGP_v4, Ensembl Plants v37，共 46,117 个基因）
- SampleID：样本编号（如：bm001 - bm467）
- ReadCount：原始读取计数
- nRC：标准化读取计数（nRC = ReadCount / sizeFactor）
- rCPM：原始 CPM（每个样本/库的总和为 1,000,000）
- rFPKM：使用 rCPM 和基因外显子长度计算的原始 FPKM
- CPM：使用 edgeR 计算的 CPM（CPM = rCPM / normFactor）
- FPKM：使用 CPM 和基因外显子长度计算的 FPKM
ta - 等位基因特异性计数表（tibble）
- SampleID：样本编号（如：bm001 - bm467）
- gid：基因 ID（AGP_v4 基因 ID）
- n0, n1：分别支持父本和母本等位基因的读取数
- ncft：支持父本和母本等位基因的冲突读取数，这种情况较少，通常是由于插入缺失区域的错配所致

biomAP 重测序数据集的短变异调用

方法

使用 [fastp] 对读取数据进行了修剪，使用 BWA 将数据比对到玉米 B73 AGP_v4 参考基因组。
对比对后的 BAM 文件进行重复标记和碱基重校准，使用 [GATK4] 工具进行处理。
针对每个样本使用 [GATK4] 的 HaplotypeCaller 调用 GVCF 文件，然后使用 [GATK4] 的 GenotypeGVCF 进行联合变异调用，生成原始的变异集。
接着，使用 [GATK4] 的 VariantRecalibration 对原始变异集进行过滤，生成最终的变异调用集。