牛轉(zhuǎn)錄組中調(diào)控變異的全面目錄

Liu, Shuli, et al. "A comprehensive catalogue of regulatory variants in the cattle transcriptome."?bioRxiv?(2020).

Abstract

表征牲畜轉(zhuǎn)錄組的遺傳調(diào)控變異體對于解釋經(jīng)濟價值特征的分子機制以及通過人工選擇提高遺傳增益的速率至關重要缓呛。在這里瞧柔，我們基于代表100多個牛組織的公開可用數(shù)據(jù)集中的11,642個RNA序列，為研究社區(qū)構(gòu)建了牛基因型組織表達圖集（cGTEx蹂匹，http：//cgtex.roslin.ed.ac.uk/）。 我們描述了跨組織轉(zhuǎn)錄組的情況戚揭，并報告了成千上萬個順式和轉(zhuǎn)基因與基因表達和24種主要組織的選擇性剪接的關聯(lián)雕旨。 我們評估了跨組織的這些遺傳調(diào)控作用的特異性/相似性，并使用多組學數(shù)據(jù)的組合在功能上對其進行了注釋悲柱。 最后锋喜，我們使用大型轉(zhuǎn)錄組關聯(lián)研究（transcriptome-wide association study，TWAS）將不同組織中的基因表達與43個重要的經(jīng)濟性狀聯(lián)系起來豌鸡，以提供新穎的生物學見識嘿般，以了解牛農(nóng)藝性狀的分子調(diào)控機制。

Introduction

①?在牛性狀相關組織中涯冠，農(nóng)藝性狀的GWAS信號表達的基因調(diào)控區(qū)域中顯著豐富炉奴，但就個體和組織數(shù)目而言，剖析基因表達變異的嘗試通常很少蛇更；

②?我們通過分析11,642種公開可用的牛RNA-Seq序列瞻赶，描述了超過100種不同的組織和細胞類型，并通過一個門戶網(wǎng)站（http://cgtex.roslin.ed.ac.uk/）使各研究部門可以自由派任，輕松地獲得結(jié)果

③?我們提出了一個標準砸逊，用于統(tǒng)一整合11,642個公共RNA-Seq數(shù)據(jù)集，并識別順式和反式表達以及剪接24個重要牛組織的定量性狀位點（eQTL和sQTL）掌逛。

--- 通過直接從RNA-Seq讀數(shù)中調(diào)用變體并使用1000 Bull Bullomes Project7的數(shù)據(jù)推算到序列水平來實現(xiàn)后者师逸，就像以前對人類數(shù)據(jù)所做的那樣;

--- 接下來，我們進行了?in silico analyses豆混，為eQTL和sQTL標注了各種公開的組學數(shù)據(jù)篓像，包括DNA甲基化，染色質(zhì)狀態(tài)和染色質(zhì)構(gòu)象特征;

--- 最后崖叫，我們通過轉(zhuǎn)錄組全關聯(lián)研究（TWAS）遗淳，將基因表達結(jié)果與來自43個牛性狀的27,214頭大型公牛的GWAS進行了整合，并鑒定了與43個性狀相關的442個先前未知的基因

Result

1心傀、Data summary

① 我們使用統(tǒng)一的質(zhì)量控制管道分析了8,642個樣品中的11,642個可公開獲得的RNA-Seq屈暗，可產(chǎn)生約2000億次清潔讀數(shù)

②?single/paired reads, clean read number, read length, sex, age, and mapping rate across samples show that the quality of publicly available data is acceptable for the analyses shown here.

③?我們還分析了21種牛組織的144種可公開獲得的全基因組亞硫酸氫鹽測序（whole-genome bisulfite sequencing , WGBS）數(shù)據(jù)集，以研究DNA甲基化的組織特異性并在功能上注釋eQTL和sQTL

Variation in gene expression across individuals and tissues

①? TPM(Transcripts Per Kilobase Million) 基因的絕對表達量，用于標準化測序深度和基因長度

②?Only 61 genes were not expressed in any of the samples. About half of those (54.10%) were located in unplaced scaffolds, with significantly (P < 0.05, 1-sided) shorter gene length, fewer exons, higher CG density, and lower sequence conservation than expressed genes

③?同樣养叛，隨著clean reads次數(shù)的增加种呐，我們檢測到了更多的 alternative splicing events

④?Furthermore, 27% of them were housekeeping RNAs and included snRNA, snoRNA, snRNAs, snoRNAs and rRNAs known to play important roles in RNA splicing

⑤?Genes without splicing events were significantly engaged in the integral component of membrane and G-protein coupled receptor signaling pathway（膜和G蛋白偶聯(lián)受體信號轉(zhuǎn)導通路）

2、Tissue specificity of gene expression

①?Tissue-specificity of gene expression was conserved across cattle and human (Fig. 2a- b), and the function of genes with tissue-specific expression accurately reflected the known biology of the tissues

②?We also calculated tissue-specificity of promoter DNA methylation and alternative splicing (Methods)

③?We found that, based on tissue-specificity, gene expression level was significantly (FDR < 0.05) and negatively correlated with DNA methylation level in promoter (Fig. 2c), and positively correlated with splicing ratios

3弃甥、Discovery of expression and splicing QTLs

①?與人類發(fā)現(xiàn)的重大變異（significant variants 爽室，eVariant）一致，這些變異以被測基因組的轉(zhuǎn)錄本起始位點（the transcript start sites淆攻，TSS）為中心阔墩。

②?我們發(fā)現(xiàn)46％（范圍14.5-73.9％）的eGenes具有一個以上與表達相關的獨立SNP（圖3c），表明基因表達的廣泛等位基因異質(zhì)性

③?等位基因特異性表達（Allele-specific expression瓶珊，ASE）分析發(fā)現(xiàn)啸箫，相關的遺傳變異在cis-eQTL中明顯過量表達，并且其效應大小顯著相關

④?Patterns and biological mechanisms underlying tissue specificity/similarity of cis-QTLs provide insights into pleiotropic regulatory effects on phenotypes（順式QTL的組織特異性/相似性背后的模式和生物學機制提供了對表型的多效調(diào)節(jié)作用的見解）

⑤?We therefore conducted a non-model-based pairwise analysis using significant eGene-eVariant pairs in one tissue to estimate the proportion of non-null associations of these pairs in another tissue (Methods).（因此伞芹，我們在一個組織中使用重要的eGene-eVariant對進行了非基于模型的成對分析忘苛，以估計在另一組織中這些對的非空關聯(lián)的比例（方法））

⑥?We speculated that the large number of trans-eQTLs in cattle could be due to the high selection intensity for economically important traits and the modest effective population size. This has led to a complex inter-chromosomal pattern of linkage disequilibrium (LD) and gene co-expression. We showed this might be the case, as we found significantly higher LD for cis-eQTL & trans-eQTL pairs than cis-eQTL & random-SNP pairs (on matched chromosomes) across all tissues (Fig. 4a).（我們推測牛中大量的反式eQTL可能是由于對重要經(jīng)濟性狀的高選擇強度和適度的有效種群規(guī)模所致。這導致了連鎖不平衡（LD）和基因共表達的復雜的染色體間模式唱较。我們證明了這可能是事實扎唾，因為我們發(fā)現(xiàn)在所有組織中，順式-eQTL和反式-eQTL對的LD顯著高于順式-eQTL和隨機-SNP對（在匹配的染色體上）（圖4a））

4南缓、Functional annotation of QTLs（QTLs的功能注釋）

① 我們采用了多層生物學數(shù)據(jù)來更好地定義遺傳調(diào)控作用的分子機制

②?如預期的那樣胸遇， cis-e/sQTLs的功能元件（例如3’UTR 和 open chromatin regions）顯著豐富

③ 與其他組織相比，橫跨13個組織的Hypomethylated regions（次甲基化區(qū)域）富含cis-e / sQTL

④?Topologically associated domains (TADs) enable chromatin interactions between distal regulatory regions and target promoters （拓撲相關結(jié)構(gòu)域（TAD）使遠端調(diào)節(jié)區(qū)和靶標啟動子之間的染色質(zhì)相互作用）

⑤?通過檢查牛肺組織的Hi-C數(shù)據(jù)汉形，我們獲得了TAD和重要的Hi-C接觸點狐榔，這些接觸點在整個組織中都是保守的15。通過與具有匹配距離的隨機eGene-SNP對進行比較获雕，我們觀察到大多數(shù)組織中TAD中eGene-eVariant對的百分比明顯更高（圖5f）。例如收捣，APCS and its cis-eQTL peak （TSS上游144 kb）被TAD包圍届案，并通過顯著的Hi-C接觸而鏈接，從而允許通過遠端（距TSS> 2 kb）調(diào)節(jié)其表達

5罢艾、eQTLs and complex trait associations

① 這項研究的主要目的是為闡明牛的農(nóng)藝性狀的遺傳和生物學基礎提供資源楣颠。因此，我們評估了在每種組織中檢測到的e / sQTL與四個不同農(nóng)藝性狀的關聯(lián)?

② TWAS增強了我們檢測因果基因并更好地了解這些特征的生物學基礎的能力

Methods

Quantification of gene expression

①?We downloaded 11,642 RNA-Seq runs (by July, 2019) from SRA (n = 11,513, https://www.ncbi.nlm.nih.gov/sra/) and BIGD databases (n = 129, https://bigd.big.ac.cn/bioproject/).

②?we first removed adaptors and low quality reads using Trimmomatic (v0.39) with parameters: adapters/TruSeq3SE.fa:2:30:10 LEADING:3 329 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36.

③?我們篩選出干凈讀數(shù)≤500K的樣品咐蚯，得到7,680個樣品童漩，并使用參數(shù)為outFilterMismatchNmax 3，outFilterMultimapNmax10和outFilterScoreMinOverLread 0.66 的STAR（v2.7.0）的單個或成對映射模塊將干凈讀數(shù)映射到ARS-UCD1.2牛參考基因組24春锋。

④?我們保留了7,264個樣本矫膨，其uniquely mapping rates≥334 60％（平均值為91.07％；范圍為60.44％-100％；表S1中的映射詳細信息）

⑤?然后侧馅，我們使用Stringtie（v2.1.1）獲得了27,608個Ensembl（v96）注釋基因的歸一化表達（TPM）危尿，并使用featureCounts（v1.5.2）提取了它們的原始讀取計數(shù)。

⑥?最后馁痴，我們使用R包dendextend中實現(xiàn)的層次聚類方法谊娇，基于log2（TPM +1）對7,264個樣本進行聚類，距離=（1-r）罗晕，其中r為Pearson相關系數(shù)济欢。我們排除了具有明顯聚類錯誤的樣本（例如，標記為肝的樣本未與其他肝樣本聚類）小渊，從而得到7,180個樣本用于后續(xù)分析法褥。

Quantification of alternative splicing

We used Leafcutter (v0.2.9)?to identify and quantify variable alternative splicing events of genes by leveraging information of junction reads (i.e., reads spanning introns) that were obtained from the STAR alignment.（有提供Leafcutter處理過程和腳本說明，此處不列舉）

Genotyping and imputation

①?我們按照基因組分析工具包（GATK）（v4.0.8.1）28中建議的最佳做法管線（默認設置）粤铭，分別將1000個Bull Bull基因組計劃中的7180個高質(zhì)量RNA-Seq樣品稱為已知基因組變異的基因型挖胃。

②最后，我們獲得了6,123個樣本梆惯，這些樣本已成功進行了基因分型和估算酱鸭。我們過濾掉了MAF <0.05和 dosage R-squared (DR2)<0.8的變體，從而產(chǎn)生了3,824,444個SNP用于QTL定位垛吗。

③?然后凹髓，我們測量了跨13個組織的WGS-SNP和RNA-Seq /估算的SNP之間的基因型一致性比率。 然后怯屉，我們使用plink（v1.90）30（--indeppairwise 1000 5 0.2）提取了153,913個 LD-independent 的SNP蔚舀，并在EIGENSOFT（v7.2.1）31中使用這些SNP對所有6,123個樣品進行了PCA分析。 我們還使用這些獨立的SNP刪除了重復的個體锨络，從而計算出樣品之間的identity- by-state distance（IBS距離）

Allele specific expression (ASE)

We conducted ASE analysis using the GATK ASEReadCounter tool (v4.0.8.1)

DNA methylation analysis of WGBS data

①?我們（2019年7月）從21種不同的組織中下載了144個可公開獲得的WGBS牛數(shù)據(jù)（表S2）赌躺。

②?我們首先使用FastQC（v0.11.2）和Trim Galore v0.4.0（--max_n 15-quality 20 -length 20 -e 0.1）分別確定讀取質(zhì)量和過濾低質(zhì)量的讀取。

③?然后羡儿，我們使用帶有默認參數(shù)的Bismark軟件（v0.14.5）32將純凈讀段映射到相同的牛參考基因組（ARS-UCD1.2）礼患。

④?After deduplication of reads, we extracted methylation levels of cytosines using the bismark_methylation_extractor (--ignore_r2 6) function.

⑤?所有WGBS數(shù)據(jù)的覆蓋率都是根據(jù)干凈的讀數(shù)計算得出的，覆蓋率為5到47倍掠归，平均覆蓋率為27.6倍缅叠。最終，我們保留了至少由5個讀數(shù)代表的CpG位點虏冻，用于后續(xù)分析肤粱。我們使用分層聚類和t-SNE方法基于共享CpG的DNA甲基化水平對樣本進行聚類

Identification of TAD and significant Hi-C contacts

①?為了發(fā)現(xiàn)遠端eVariant與目標eGenes之間的潛在染色質(zhì)相互作用，我們從牛的肺組織的Hi-C數(shù)據(jù)中鑒定了TAD和Hi-C接觸物

②?我們使用Trim Galore（v0.4.0）修整了adapter sequences 和低質(zhì)量讀值（--max_n 15 --quality 20 --length 20 -e 0.1）厨相，產(chǎn)生了約8.2億次純凈讀取领曼。

③?然后鸥鹉，我們使用BWA將純凈讀圖映射到牛參考基因組（ARS-UCD1.224）。我們應用HiCExplorer v3.4.1構(gòu)建了一個分辨率為10kb的 Hi-C contacts 矩陣悯森，并使用hicFindTAD識別了TAD宋舷。

④?我們使FDR小于0.01的TAD保持eQTL與eGenes的鏈接。我們進一步采用HiC-Pro（v2.11.4）從Hi-C數(shù)據(jù)中以10 kb分辨率調(diào)用Hi-C contacts

⑤?簡而言之瓢姻，HiC-Pro使用Bowtie2（v2.3.5）34對牛參考基因組進行了比對的純凈讀取祝蝠。建立a contact matrix后，HiC-Pro生成染色體內(nèi)和染色體間圖譜幻碱，并使用原始的ICE歸一化算法對其進行歸一化绎狭。我們認為FDR <0.05的Hi-C接觸很重要。

Tissue-specificity analysis of gene expression, alternative splicing and DNA methylation

①?為了量化基因的組織特異性表達褥傍，我們計算了114個組織中每個基因的 t-統(tǒng)計量儡嘶。

②?We scaled the log2-transformed expression (i.e., log2TPM) of genes to have a mean of zero and variance of one within each tissue.

③?然后，我們使用最小二乘法（least-squares）擬合每個組織中每個基因的線性模型

④?為了檢測組織特異性的選擇性剪接恍风，我們使用 leafcutter 通過比較目標組織和其余組織的樣本來分析差異內(nèi)含子切除

⑤?我們使用Benjamini-Hochberg方法（FDR）來控制多重測試

⑥ 對于DNA甲基化蹦狂，我們著眼于基因啟動子（從TSS的上游1500bp到下游500bp）和身體區(qū)域（從TSS到TES）的DNA甲基化水平，這是使用roimethstat函數(shù)MethPipe（v3.4.3）通過加權(quán)甲基化方法計算的

⑦ 我們使用與組織特異性表達分析相同的方法朋贬，計算了每個基因的啟動子的 t-統(tǒng)計量凯楔。

⑧?我們還使用參數(shù)-t DeNovoDMR -MR 0.5 -AG 1.0 -MS 0.5 -ED 0.2 -SM 0.6-CD 500 -CN 5-SL 20 -PD 0.05-使用SMART2在全基因組模式下檢測了組織特異性甲基化區(qū)域

Covariate analysis for QTL discovery

為了解決基因表達中轉(zhuǎn)錄組范圍內(nèi)變異( transcriptome-wide variation )的 hidden batch effects 和其他技術/生物學來源，我們使用表達殘基概率估計（PEER）方法估算了每個組織中的潛在協(xié)變量锦募。

cis-eQTL mapping

Meta-analysis of cis-eQTLs of muscle samples from three sub-species

We then conducted a meta-analysis to integrate cis-eQTL results from three sub-species using the Metal tool45.

cis-sQTL mapping

在這24個組織的每一個中摆屯，我們應用了在FastQTL41中實施的線性回歸模型，以測試目標內(nèi)含子簇上下游1 Mb內(nèi)的基因型（MAF> 1％）及其對應的內(nèi)含子切除率的關聯(lián)糠亩。

trans-eQTL mapping

WGCNA co-expression network analysis and estimation of π1 statistics

We applied the Weighted Gene Co-Expression network (WGCNA) analysis to obtain gene co-expression network in each of 24 tissues used in eQTL mapping. We estimated the sharing of cis-eQTLs and cis-sQTLs among tissues using the π1 statistics, as described in human GTEx Consortium (2015).

TWAS analyses

①?To associate gene expression in a tissue with complex traits, we conduced TWAS analysis using S-PrediXcan48 by prioritizing GWAS summary statistics for 43 agronomic traits of economic importance in cattle, including reproduction (n = 11), production (milk- relevant; n = 6), body type (n = 17), and health (immune/metabolic-relevant; n = 9).（為了使基因表達與具有復雜性狀的組織相關聯(lián)虐骑，我們通過優(yōu)先利用GWAS摘要統(tǒng)計數(shù)據(jù)優(yōu)先考慮牛的43個具有重要經(jīng)濟價值的農(nóng)藝性狀，包括繁殖（n = 11）赎线，生產(chǎn)（與牛奶相關廷没； n = 6），體型（n = 17）和健康（免疫/代謝相關垂寥； n = 9）腕柜。 ）

②?For body conformation (type), reproduction, and production traits, we conducted a single-marker GWAS by fitting a linear mixed model in 27,214 U.S. Holstein bulls as described previously.(對于身體構(gòu)象（類型），繁殖和生產(chǎn)特征矫废，我們通過擬合線性混合模型在27,214個美國荷斯坦公牛中進行了單標記GWAS，如前所述)

③?我們使用基因型和表達數(shù)據(jù)構(gòu)建了一個嵌套的交叉驗證的彈性網(wǎng)預測模型(a Nested Cross Validated Elastic Net prediction model)砰蠢。

④?We visualized the Manhattan plots of P-values of all tested genes using ggplot2 (v3.3.2) in R（我們使用R中的ggplot2（v3.3.2）可視化了所有測試基因的P值的曼哈頓圖 ）

Other downstream bioinformatic analysis（其他下游生物信息學分析）

①?We used Genomic Association Tester (GATv1.3.4) 10,000 permutations to estimate the functional enrichment of QTLs in particular genomic regions, chromatin states, methylation elements, and WGCNA co-expression modules.（我們使用了基因組關聯(lián)測試儀（GATv1.3.4）10,000個排列來估計特定基因組區(qū)域蓖扑，染色質(zhì)狀態(tài)，甲基化元件和WGCNA共表達模塊中QTL的功能富集台舱。） We considered enrichments with FDR (adjusted P-values with Benjamini-Hochberg method) < 0.05 as significant.

② We used the R package, ClusterProfiler, to annotate function of genes based on the Gene ontology database from Bioconductor (org.Bt.eg.db v3.11.4). We considered GO terms with FDR < 0.05 as significant. （我們使用R包ClusterProfiler律杠，基于Bioconductor（org.Bt.eg.db v3.11.4）的基因本體數(shù)據(jù)庫對基因的功能進行注釋潭流。我們認為FDR <0.05的GO項很重要。）

③?We obtained complex traits/diseases that were associated with a particular gene in humans using the Region PheWAS function in GeneATLAS database (http://geneatlas.roslin.ed.ac.uk/) with set region of ±50kb and P-value threshold of 10-8 （我們使用GeneATLAS數(shù)據(jù)庫（http://geneatlas.roslin.ed.ac.uk/）中的Region PheWAS函數(shù)獲得了與人類特定基因相關的復雜性狀/疾病柜去，其設定區(qū)域為±50kb灰嫉，P值閾值為 的10-8）

Data?availability statement?

本研究中分析的所有原始測序數(shù)據(jù)均可在NCBI Gene Expression Omnibus（GEO;?https://www.ncbi.nlm.nih.gov/geo/）中公開獲得。所有已處理的結(jié)果和腳本代碼都可以在?https://cgtex.roslin.ed.ac.uk/?中找到嗓奢。