Transcriptomic characteristics of bronchoalveolar lavage fluid and peripheral blood mononuclear cells in COVID-19 patients
文章： 新冠肺炎患者支氣管肺泡灌洗液與外周血單核細胞的轉(zhuǎn)錄組學特征
期刊： SSRN
影響因子： NA
材料： COVID-19患者的支氣管肺泡灌洗液和外周血單核細胞分離的RNA
技術(shù)： RNA-seq
分享人： 賴為潔
摘要
COVID-19持續(xù)爆發(fā)，對公共衛(wèi)生構(gòu)成了巨大威脅。但是，我們對這種致命疾病的宿主免疫反應的了解有限，是我們無法判斷有效的支持治療藥物和治療方法。
為了獲取宿主對SARS-CoV-2（HCoV-19）感染的炎癥反應的轉(zhuǎn)錄特征，我們對從COVID-19患者的支氣管肺泡灌洗液（BALF）和外周血單核細胞（PBMC）標本中分離的RNA進行了轉(zhuǎn)錄組測序。
研究結(jié)果揭示了患者SARS-CoV-2感染的不同宿主炎癥細胞因子譜，并揭示了COVID-19發(fā)病與細胞因子過度釋放（例如CCL2 / MCP-1，CXCL10 / IP-10，CCL3 / MIP-1A以及CCL4 / MIP1B）之間的關(guān)聯(lián)。
此外，SARS-CoV-2誘導淋巴細胞凋亡和P53信號通路的激活可能是患者淋巴細胞減少的原因。
COVID-19患者的轉(zhuǎn)錄組數(shù)據(jù)集將為抗炎藥物的臨床指導和了解宿主反應的分子機制提供寶貴的資源。
材料
3例SARS-CoV-2患者和3例健康供體獲得了全血和PBMC（Ficoll制劑）
方法
1.RNA-seq reads were firstly mapped to rRNA sequences to remove potential rRNA reads using STAR (v2.7.2b)with the default parameter.
首先使用默認參數(shù)STAR（v2.7.2b）將RNA-seq reads 映射到rRNA序列，以去除潛在的rRNA reads。（寡聚-dT（oligo-dTs ）共價偶聯(lián)的磁珠純化信使RNA）
2.The rest reads were then mapped to the human genome (hg38) with GENCODE gene annotation (v32) with parameters “–sjdbScore 1 –outFilterMultimapNmax?20 –outFilterMismatchNmax?999 –outFilterMismatchNoverReadLmax?0.04 –alignIntronMin 20 –alignIntronMax?1000000 –alignMatesGapMax?1000000 –alignSJoverhangMin 8 –alignSJDBoverhangMin 1” following the guideline of ENCODE RNA-seq pipeline (https://github.com/ENCODE-DCC/long-rna-seq-pipeline).
將其余的reads通過GENCODE基因注釋（v32）映射到人類基因組（hg38），參數(shù)為“ –sjdbScore 1 –outFilterMultimapNmax 20 –outFilterMismatchNmax 999 –outFilterMismatchNoverReadLmax 0.04 –alignIntronMin 20 –alignIntronMax 1000000 –alignMatesGapMax 1000000 –alignSJoverhangMin 8 ”遵循ENCODE RNA-seq管道指南（https://github.com/ENCODE-DCC/long-rna-seq-pipeline）。
- The human un-mappable reads were then mapped to the SARS-CoV-2 genome as previously reported. PCR replicates mapped in the human genome were removed with picard MarkDuplicates program (v2.13.2-1) . RNA-seq signal tracks were generated by using bam2wig.py provided in RSeQC package (v3.0.0), and visualized in UCSC genome browser with custom track hub.
然后，將不能比對到人的reads定位到SARS-CoV-2基因組。用picard MarkDuplicates程序（v2.13.2-1）刪除了映射在人類基因組中的PCR復制。
通過使用RSeQC軟件包（v3.0.0）中提供的bam2wig.py生成RNA-seq信號軌道，并在UCSC基因組瀏覽器中使用自定義軌道中心對其進行可視化。
4.Gene expression was calculated by featureCounts in SubReads package (v1.5.3) with the “-M” parameter. Differentially expressed genes were called by using DESeq2 package (v1.26.0) with the following criteria: adjusted p-value?<?0.05 and fold-change?>?2. The expressed genes (requiring reads counts greater than 10 and 100 for BALF and PBMC, respectively), were selected and normalized to counts per millions (CPM) for further analysis.
基因表達是通過SubReads軟件包（v1.5.3）中的featureCounts 使用“ -M”參數(shù)來計算的。
通過使用DESeq2程序包（v1.26.0）調(diào)用差異表達的基因，其標準如下：調(diào)整后的p值<0.05，倍數(shù)變化>2。
選擇表達的基因（分別要求BALF和PBMC的讀數(shù)分別大于10和100），并將其標準化為每百萬計數(shù)（CPM）進行進一步分析。
5.Functional enrichment analysis was performed on the list of differentially expressed genes by using the clusterProfiler package (v 3.14.3)to determine if the genes are enriched for specific terms. The level of significance for the enrichment was calculated by a hypergeometric test for each term using all expressed genes as background. The p-values were corrected for multiple hypothesis tests using the Benjamini and Hochberg method to control the false discovery rate (FDR). To inspect genes in specific GO-term and pathway, the GO and KEGG pathway annotation were downloaded from Gene Ontology Resource (Last updated: 22 February 2020) and KEGG database (Last updated: 14 January 2020) , respectively. The interacting proteins to the differentially expressed genes in PBMC or BALF were identified from the functional human network InWeb_InBipMap (a human protein–protein interaction network with direct interaction evidence in high-confidence). Cytoscape (3.5.0) was used for visualizing the PPI subnetwork.
通過使用clusterProfiler程序包（v 3.14.3） 對差異表達的基因列表進行功能富集分析，以確定基因是否針對特定術(shù)語富集。
通過超幾何學測試，使用所有表達的基因作為背景，計算每個項目的富集顯著性水平。
使用Benjamini和Hochberg方法控制錯誤發(fā)現(xiàn)率（FDR），對多個假設檢驗的p值進行了校正。要檢查特定GO術(shù)語和途徑中的基因，可從Gene Ontology資源（最新更新：2020年2月22日）和KEGG數(shù)據(jù)庫（最新更新：2020年1月14日 下載GO和KEGG途徑注釋。
PBMC或BALF中與差異表達基因的相互作用蛋白是從功能性的人為網(wǎng)絡InWeb_InBipMap（具有高可信度的直接相互作用證據(jù)的人蛋白質(zhì)-蛋白質(zhì)相互作用網(wǎng)絡）中鑒定出來的。Cytoscape（3.5.0）用于可視化PPI子網(wǎng)。
可從https://github.com/zhouyulab/ncov/獲得分析的源代碼
結(jié)果:
我們對RNA-seq數(shù)據(jù)進行了質(zhì)量控制分析，并分別在補充表1和2的BALF和PBMC樣本中總結(jié)了映射到人類基因組和SARS-CoV-2基因組（GenBank MN988668）的讀數(shù)的統(tǒng)計信息。
與我們先前針對病毒基因組組裝的報告[ 3 ]一致，BALF患者樣品包含高比例的病毒讀數(shù)，而健康對照數(shù)據(jù)集的計數(shù)為零。
有趣的是，來自患者的PBMC樣本幾乎沒有顯示出病毒讀數(shù)，表明SARS-CoV-2可能不會感染PBMC。數(shù)據(jù)說明了對照組或患者組之間的高度一致性，分別通過BALF和PBMC樣本的相關(guān)性和聚類分析得到證明（補充圖S1A和S1B）。
有趣的是，如基因組瀏覽器中顯示的RNA-seq信號所示，三分之二的患者IL6和IL6R降低（D，E），盡管TP53表達顯示出增加的趨勢（F），盡管它們在統(tǒng)計學上沒有被確定為顯著改變的基因。
調(diào)控基因的功能富集分析
作者對支氣管肺泡灌洗液和外周血單核細胞中差異調(diào)節(jié)基因進行了基因功能富集分析。在支氣管肺泡灌洗液樣本中，上調(diào)的基因與病毒的入侵有關(guān)，病毒感染引起各種膜結(jié)構(gòu)和內(nèi)質(zhì)網(wǎng)的變化；而在外周血單核細胞樣本中上調(diào)的基因主要富含“補體激活”、“循環(huán)免疫球蛋白介導的體液免疫反應和“ B細胞介導的免疫”，表明外周血單核細胞中的免疫活性已激活；此外，還激活了一系列與炎癥相關(guān)的過程。
進一步對這些上調(diào)和下調(diào)的基因進行了KEGG通路分析。在支氣管肺泡灌洗液樣本中，富集改變基因的途徑包括'核糖體'，'內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工'，'吞噬體'，'戊糖磷酸途徑'，'碳代謝'和'溶酶體'。對于外周血單核細胞樣本，“細胞周期”和“ p53信號傳導”途徑豐富了上調(diào)的基因，而與細胞因子相關(guān)的途徑則豐富了下調(diào)的基因。
調(diào)節(jié)細胞因子的分析
COVID-19患者在SARS-CoV感染會發(fā)生細胞因子風暴。我們鑒定了BALF和PBMC中所有表達的細胞因子，其中帶有星號標記的基因發(fā)生了顯著變化。如在臨床病例中，在BALF中觀察到高細胞因子表達。與實驗室檢查結(jié)果一致，我們發(fā)現(xiàn)細胞因子IL10，CCL2 / MCP-1，CXCL10 / IP-10，CCL3 / MIP-1A和CCL4 / MIP1B在患者的BALF樣品中高度表達。
兩名患者（WHU01和WHU02）之間的進一步比較表明，某些細胞因子的表達在患者之間是可變的，例如IL10，IL36RN，IL36G，TNFSF15，CCL5，TNFSF10，CXCL1和IL33。
SARS-CoV-2感染可能導致淋巴細胞凋亡
3名患者的實驗室檢查結(jié)果（補充表3）表明，包括患者血液中的淋巴細胞在內(nèi)的各種類型免疫細胞的細胞計數(shù)減少。先前的臨床和尸檢報告表明，患者的淋巴細胞已大大減少。
我們發(fā)現(xiàn)一些顯著改變的基因豐富了細胞凋亡和P53信號通路（[圖5]，包括CTSL，CTSB，DDIT4，RRAS，CTSD，BIRC5，TNFSF10，CTSZ，NTRK1，IGFBP3，CCNB1，RRM2，CCNB2，GTSE1， CDK1，STEAP3和TP53I3。有趣的是，TP53是凋亡過程中的重要基因，在兩名患者中顯示出增加的趨勢，這表明PMBC的減少可能是由于凋亡引起的。
比較患者BALF和PBMC中不同變化的基因
基于差異基因表達分析，我們發(fā)現(xiàn)BALF樣品和PBMC之間的細胞內(nèi)基因變化水平存在顯著差異。這種區(qū)別是由于對兩種類型細胞的病毒感染不同。
值得注意的是，該病毒似乎并未感染PBMC，因為我們未能檢測到PBMC中的病毒RNA和ACE2表達，這也得到了骨髓血細胞表達數(shù)據(jù)庫的支持。
為了說明兩種細胞類型之間的差異，我們對BALF樣品和PBMC之間具有不同趨勢的基因進行了分類（A，B））。
有17個基因朝同一方向變化，而36個基因則有相反的趨勢。
我們首先分析了PBMC中增加但BLAF中減少的基因。在這15個基因中，有5個基因參與了小分子分解代謝過程：ADA2，HK1，MGAT1，PGD和PLA2G15。中性粒細胞的激活和免疫涉及4個基因：ADA2，CTSD，GAA，LAIR1。接下來，我們分析了在BLAF中增加但在PBMC中減少的基因，這些基因在[圖6]（B）中列出。盡管兩個樣品之間有17個基因以相同的趨勢變化，但生物學意義仍需進一步研究。
