基于Linux平臺構(gòu)建單細(xì)胞轉(zhuǎn)錄組病毒信息分析系統(tǒng)*

陸庭希，侯晨陽，謝 興，胡艷玲

（廣西醫(yī)科大學(xué)生命科學(xué)研究院，南寧 530021）

隨著測序技術(shù)的飛速發(fā)展，分析處理測序數(shù)據(jù)信息已然成為現(xiàn)代生物研究的常規(guī)手段，然而傳統(tǒng)分析流程過于碎片化，步驟繁瑣且耗時較長。生物信息分析系統(tǒng)包括生物數(shù)據(jù)預(yù)處理和分析工具兩大宏觀部分[1]。通常所指的生物信息分析系統(tǒng)是包括收集并進(jìn)行預(yù)處理相關(guān)生物數(shù)據(jù)（包括核酸序列、蛋白質(zhì)序列、蛋白質(zhì)結(jié)構(gòu)或知識庫），集成自主開發(fā)軟件、腳本以及相關(guān)的公共開源的分析工具，可供用戶進(jìn)行數(shù)據(jù)處理、數(shù)據(jù)分析、生物信息挖掘等功能的統(tǒng)一分析平臺[2]。針對當(dāng)前數(shù)據(jù)量大、關(guān)聯(lián)復(fù)雜、方法多樣、工具軟件繁多等突出的需求現(xiàn)狀，構(gòu)建整合的數(shù)據(jù)分析系統(tǒng)，集成已有的平臺及特定軟件集合，達(dá)到加速分析速度、減少分析過程中隨機(jī)誤差，降低分析流程復(fù)雜度的效果，顯得尤為重要。現(xiàn)階段常規(guī)分析方法仍以各個軟件單獨(dú)調(diào)參為主要手段，分析流程中包含的所有軟件均需要使用者手工修改相關(guān)代碼參數(shù)，保證其正確有效運(yùn)行來得出期望的結(jié)果。局限性在于目前的分析方法代碼復(fù)現(xiàn)率低，在不同的設(shè)備、不同的系統(tǒng)環(huán)境下、不同的數(shù)據(jù)集中難以復(fù)現(xiàn)，對于大規(guī)模、多樣本數(shù)據(jù)分析，存在分析效率低、隨機(jī)誤差大、分析結(jié)果不穩(wěn)定等問題。生物信息分析流程自動化、標(biāo)準(zhǔn)化是亟待解決的問題。因此，基于Linux 平臺與單細(xì)胞轉(zhuǎn)錄組測序數(shù)據(jù)，整合現(xiàn)階段分析流程所使用的相關(guān)分析軟件，將其無縫串聯(lián)起來形成一個標(biāo)準(zhǔn)化的快速分析流程，建立該病毒分析信息系統(tǒng)是非常必要的，以助于達(dá)到快速分析醫(yī)學(xué)大數(shù)據(jù)的目的。本文以Linux 系統(tǒng)為基礎(chǔ)，整合多個生物信息分析軟件，輔以自動化腳本和相關(guān)數(shù)據(jù)庫構(gòu)建單細(xì)胞轉(zhuǎn)錄組病毒信息分析系統(tǒng)，現(xiàn)將結(jié)果報(bào)道如下。

1 單細(xì)胞轉(zhuǎn)錄組病毒信息分析系統(tǒng)構(gòu)建過程

1.1 分析系統(tǒng)環(huán)境配置

本系統(tǒng)基于Linux（Version:3.10.0-693.el7.x86_64,builder@kbuilder.dev.centos.org）操作系統(tǒng)開發(fā)，硬件使用多臺超級計(jì)算機(jī)（CPU:32 個Intel(R)Xeon(R)Silver 4110 CPU@2.10GHz）組成的分布式服務(wù)器/物理計(jì)算節(jié)點(diǎn)使用DDR4 2666MHz內(nèi)存125G或250G/70TB 冗余硬盤陣列（RAID 1 方式），軟件部分使用開源的Linux 系統(tǒng)Centos，并進(jìn)行基本的系統(tǒng)設(shè)置及環(huán)境配置，分析軟件的安裝及環(huán)境配置。并且，由于科學(xué)計(jì)算任務(wù)對系統(tǒng)的穩(wěn)定性和可靠性提出了較高的要求，所以分布式服務(wù)器安裝在室溫穩(wěn)定且低于25oC 的機(jī)房環(huán)境下，保證其正常工作和運(yùn)行[3]。環(huán)境配置如下：export PATH=$PATH:/public/software/sratoolkit/bin/;export PATH=$PATH:/$User_PATH/prepare_fastq_v1/;export PATH=$PATH:/public/software/samtools/bin/;source/public/software/ellranger-3.1.0/sourceme.bash

1.2 系統(tǒng)分層設(shè)計(jì)流程

單細(xì)胞轉(zhuǎn)錄組病毒信息分析系統(tǒng)采用分布式分層設(shè)計(jì)，其中包括上游分析層，下游分析層。每層內(nèi)整合了相關(guān)生物信息分析軟件及自動化腳本，見圖1。其優(yōu)勢在于，各層相互獨(dú)立，在分布式集群中每個模塊可以在不同計(jì)算節(jié)點(diǎn)中處理任務(wù)，便于拓展以及防止過多任務(wù)投入導(dǎo)致的宕機(jī)。同時，這種分層設(shè)計(jì)不僅僅在于物理層面的劃分，而包含邏輯層面的劃分。不同層的模塊可以在同一物理計(jì)算節(jié)點(diǎn)中運(yùn)行，而同層模塊中運(yùn)行也可以由不同的物理計(jì)算節(jié)點(diǎn)來共同并行完成[4]。

圖1 單細(xì)胞轉(zhuǎn)錄組病毒信息分析系統(tǒng)分層流程圖

1.2.1 上游分析層 結(jié)合系統(tǒng)分層流程圖，對系統(tǒng)中的各模塊進(jìn)行介紹，各個模塊集成軟件，見表1。

數(shù)據(jù)源模塊：系統(tǒng)數(shù)據(jù)來源，單細(xì)胞轉(zhuǎn)錄組測序數(shù)據(jù)來源可以是原始下機(jī)數(shù)據(jù)，同時支持GEO數(shù)據(jù)庫（gene expression omnibus database）中的公共測序數(shù)據(jù)，數(shù)據(jù)通過SFTP 文件服務(wù)器進(jìn)行上傳及存儲。該模塊整合NCBI 中的sratoolkit:fastq-dump[5]和Bamtofastq-1.2.0，進(jìn)行數(shù)據(jù)預(yù)處理。

數(shù)據(jù)處理模塊：此模塊應(yīng)完成對源數(shù)據(jù)的上游分析處理，符合預(yù)處理要求的測序數(shù)據(jù)傳入此模塊。此模塊主要集成了Cellranger、BBMap、Krak‐en2、Bracken、KrakenTools、BLAST[6]以及自主開發(fā)編寫的Python腳本kraken2barcodem.py。

1.2.2 下游分析層 根據(jù)上游分析層獲得的ma‐trix、features、barcodes 以及所提取對應(yīng)物種的bar‐code序列，通過R語言的分析，可以確定比對所得的物種所在的細(xì)胞種類。該層集成了相關(guān)的R 包，來構(gòu)建下游分析功能。R 包：Seurat（https://satijalab.org/seurat/）[7]、ggplot2、cowplot、Matrix、dplyr。

數(shù)據(jù)讀入模塊：此模塊加載集成分析所需R包，并使用CreateSeuratObject 函數(shù)創(chuàng)建Seurat 對象，并提示用戶對project/min.cell/min.featuresd 參數(shù)進(jìn)行設(shè)置。數(shù)據(jù)讀入支持10X 單細(xì)胞測序數(shù)據(jù)格式和原始表達(dá)矩陣格式的數(shù)據(jù)。

數(shù)據(jù)質(zhì)控模塊：數(shù)據(jù)質(zhì)控是至關(guān)重要的環(huán)節(jié)，直接影響后續(xù)分析數(shù)據(jù)的可靠性。Seurat質(zhì)控統(tǒng)計(jì)信息中包含每個細(xì)胞中的UMI數(shù)量和基因數(shù)，以及每個細(xì)胞中的線粒體基因占比；統(tǒng)計(jì)基因數(shù)、RNA、線粒體基因分布，并計(jì)算基因數(shù)與線粒體基因以及RNA 數(shù)量的分布相關(guān)性；篩選高質(zhì)量細(xì)胞，盡量保證線粒體基因占比低于10%，表達(dá)的基因數(shù)至少大于200，測到的RNA 分子數(shù)在1 000～20 000 之間。以上質(zhì)控結(jié)果以表圖形式存儲。

數(shù)據(jù)分析模塊：此模塊完成表達(dá)量的標(biāo)準(zhǔn)化、歸一化與PCA、確定細(xì)胞類群分析PC、細(xì)胞分群、計(jì)算marker 基因的分析。表達(dá)量的標(biāo)準(zhǔn)化使用Log‐Normalize 算法，F(xiàn)indVariableFeatures 計(jì)算表達(dá)量變化顯著的基因；歸一化消除特定變量對整體數(shù)據(jù)的影響，如批次效應(yīng)，線粒體數(shù)量的差異等。

表1 分析軟件、R包的來源與功能

1.3 整合軟件功能具體描述

上游分析部分：Sratoolkit由NCBI 提供，用于處理來自SRA 數(shù)據(jù)庫中高通量測序數(shù)據(jù)的工具包，prefetch 模塊可以批量下載數(shù)據(jù)庫中數(shù)據(jù)集的原始數(shù)據(jù)；fastq-dump 模塊可以將二代測序sra 格式的數(shù)據(jù)轉(zhuǎn)換為分析所需的fastq 格式，其中雙端測序文件使用--split-3 參數(shù)；Bamtofastq 是samtools 軟件中的功能模塊，其功能是從序列比對結(jié)果BAM 文件轉(zhuǎn)換為fastq格式文件，對于雙端測序數(shù)據(jù)需要使用-fq與-fq2 參數(shù)；Cellranger 是10X genomics 開發(fā)專門用于單細(xì)胞轉(zhuǎn)錄組測序分析的數(shù)據(jù)分析軟件，可以完成數(shù)據(jù)拆分（cellranger mkfastq）、細(xì)胞定量（cell‐ranger count）、組合分析（cellranger aggr）、參數(shù)調(diào)整（cellranger reanalyze）等相關(guān)功能，將測序數(shù)據(jù)比對至參考基因組，完成細(xì)胞和基因表達(dá)定量并生成細(xì)胞-基因表達(dá)矩陣，基于此矩陣進(jìn)行細(xì)胞聚類和差異表達(dá)分析；BBMap 是一種用于DNA 和RNA 測序reads 的比對工具，可以完成感知全局的拼接工作，它對基因組大小和contigs 數(shù)量沒有限制，并且支持多線程來達(dá)到極快的索引速度；Kraken2 是一個基于k-mer 精確匹配算法的高精度序列分類軟件，能夠快速地將測序reads進(jìn)行物種分類，進(jìn)而注釋出分析得到相關(guān)物種信息；Bracken（Bayesian Reestima‐tion of Abundance with KrakEN）是一種高精度統(tǒng)計(jì)方法，結(jié)合Kraken2 可以實(shí)現(xiàn)高準(zhǔn)確度的測序數(shù)據(jù)物種分類分析，并計(jì)算物種豐度；KrakenTools 是一套適用于Kraken2/Bracken結(jié)果后續(xù)分析的腳本，在系統(tǒng)中整合了從kraken 結(jié)果中提取指定序列的功能；BLAST（Basic Local Alignment Search Tool）是在蛋白質(zhì)數(shù)據(jù)庫或DNA 數(shù)據(jù)庫中進(jìn)行相似性比較的分析工具，BLAST能迅速與公開數(shù)據(jù)庫進(jìn)行相似性序列比較，從而得出Coverage 值和Identical 值，BLAST結(jié)果中的得分是對相似性的統(tǒng)計(jì)說明，可以完成指定序列的特異性檢驗(yàn)。

下游分析部分：Seurat 是一個R 包，其功能為識別和解釋單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)中的異質(zhì)性來源、整合了不同類型的單細(xì)胞數(shù)據(jù)的函數(shù)、單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)中的細(xì)胞質(zhì)控、聚類、差異分析、細(xì)胞注釋、以及結(jié)果可視化；ggplot2 是用于繪圖的R 拓展包，整合在系統(tǒng)中為下游分析結(jié)果提供數(shù)據(jù)可視化；cowplot是一個繪圖插件，可以為ggplot2提供出版級別的主題，優(yōu)化可視化結(jié)果；dplyr 是一個用于數(shù)據(jù)清洗和處理的R包，整合到系統(tǒng)中，為輸入數(shù)據(jù)的分析處理提供支持。

2 系統(tǒng)功能實(shí)現(xiàn)及測試結(jié)果分析

2.1 數(shù)據(jù)收集與預(yù)處理

當(dāng)前信息分析系統(tǒng)整合了數(shù)據(jù)收集與預(yù)處理所需的軟件，部署于服務(wù)器中，對軟件完成封裝。將sra、bam 格式測序數(shù)據(jù)轉(zhuǎn)換為fastq 格式數(shù)據(jù)，并使用prepare_fastq.py 提取sra 格式數(shù)據(jù)的前100 行、bam 格式數(shù)據(jù)的前1000 行進(jìn)行測試，判斷其可用性，當(dāng)數(shù)據(jù)不可用時，彈出提示并中斷程序。對于sra 格式數(shù)據(jù)：如果spot 數(shù)為1，則判斷該數(shù)據(jù)不可用；如果spot數(shù)為2，則spot1為read1，spot2為read2；如果spot 數(shù) 為3，則spot1 為index，spot2 為read1，spot3 為read2。同時檢查read1、read2 長度，如果長度低于16，則不是barcode，判斷該數(shù)據(jù)不可用；read2 長度要求不低于50 bp。對于bam 格式數(shù)據(jù)：檢查tags CB，CR 是否完整，如果CR，CB 低于一定比例，則判斷該數(shù)據(jù)不可用。符合預(yù)處理標(biāo)準(zhǔn)的數(shù)據(jù)，方能傳入下一模塊進(jìn)行分析。

2.2 分析系統(tǒng)功能分析

2.2.1 參考庫下載與用戶數(shù)據(jù)管理 在構(gòu)建此分析系統(tǒng)時，需要高效整合數(shù)據(jù)分析所用的參考庫與參考基因組，以達(dá)到系統(tǒng)高效分析的效果。使用Cellranger 所需的人類參考基因組GRCh38，從Cell‐ranger 官 網(wǎng)（https://cf.10xgenomics.com/supp/cellexp/refdata-gex-GRCh38-2020-A.tar.gz）上獲取已構(gòu)建好的索引文件；BBMap 所使用的參考基因組為Homo_sapiens.GRCh38.dna.prmary_assembly.fa.gz（ftp://ftp.ensembl.org/pub/release-101/fasta/homo_sa‐piens/dna/）；Kraken2 所使用的參考庫為minikrak‐en2_v1_8 GB（http://ccb.jhu.edu/software/kraken2/in‐dex.shtml?t=downloads）。用戶可依據(jù)不同任務(wù)自建工作文件夾在系統(tǒng)目錄中，分析結(jié)果默認(rèn)保存在系統(tǒng)目錄中，若用戶自建目錄，則存放于用戶指定目錄下。

2.2.2 系統(tǒng)數(shù)據(jù)分析 在構(gòu)建的單細(xì)胞轉(zhuǎn)錄組病毒分析系統(tǒng)中主要完成兩個主要任務(wù)：一是從GEO數(shù)據(jù)庫中獲取單細(xì)胞測序數(shù)據(jù)，轉(zhuǎn)換格式為fastq，完成質(zhì)控等數(shù)據(jù)預(yù)處理，并與參考基因組比對去除宿主，在此基礎(chǔ)上，聯(lián)合Kraken2 與Bracken 對測序reads 進(jìn)行物種注釋，尋找病毒reads，并提取指定病毒序列和barcode序列，進(jìn)行Blast驗(yàn)證其特異性，完成上游分析；二是對所獲得的matrix、features、bar‐codes 以及所提取的對應(yīng)物種的barcode 序列，通過R 語言分析來確定比對的物種所在的細(xì)胞種類，對細(xì)胞進(jìn)行分群并計(jì)算marker 基因。其中，圖2 展示了標(biāo)準(zhǔn)化之后的整體表達(dá)水平。圖3展示了PCA分析從大量的基因表達(dá)信息中，提取了對整體基因表達(dá)量影響最大的效應(yīng)，命名為PC1，PC2，基因之間的表達(dá)差異表現(xiàn)在PC1、PC2數(shù)值上的差異。圖4以2D 點(diǎn)圖展示了PCA 的結(jié)果，其中每個點(diǎn)代表了每個細(xì)胞。圖5 用熱圖展示了PCA 的結(jié)果。圖6 為聚類分析后的TSNE 的分群展示，其中每個顏色代表了一個cluster 后鑒定到的一種細(xì)胞群，散點(diǎn)代表了每個細(xì)胞，圖中數(shù)字則代表了該群的cluster編號，由圖中可知共分為了8 個不同的細(xì)胞群。圖7 展示了子群中的marker 基因的表達(dá)水平，能更為直觀地看到marker 基因在細(xì)胞中的表達(dá)水平，其中顏色越深則表示該基因在這些子群或細(xì)胞中的表達(dá)越高。圖8使用小提琴圖展示這幾個marker基因在子群中的相對表達(dá)水平，可以衡量該基因作為子群marker基因的特異性。圖9以熱圖的形式展示了marker基因在所有子群中的表達(dá)水平，并按照cluster 編號進(jìn)行排序，圖中的顏色表示了表達(dá)水平，黃色表示高表達(dá)。

圖2 表達(dá)量標(biāo)準(zhǔn)化

圖3 PCA結(jié)果-1

圖4 PCA結(jié)果-2

圖5 PCA結(jié)果-3

圖6 TSNE細(xì)胞分群

圖8 單個標(biāo)記基因子群表達(dá)水平小提琴圖

圖9 Top marker基因熱圖

3 討論

本文整合了多個可免費(fèi)獲得的生物信息學(xué)分析軟件，輔以自主開發(fā)的自動化腳本，在Linux 平臺上構(gòu)建了可完成大數(shù)據(jù)、多任務(wù)單細(xì)胞轉(zhuǎn)錄組病毒信息分析的系統(tǒng)。應(yīng)用該系統(tǒng)，我們已成功分析并獲取到多組測序數(shù)據(jù)中病毒信息，得到了許多有意義的結(jié)果。

通過對本次構(gòu)建的單細(xì)胞轉(zhuǎn)錄組病毒信息分析系統(tǒng)的深入研發(fā)，在不同生信分析軟件之間建立了接口和構(gòu)建了流程，給用戶提供了完整、易用的分析平臺，幫助用戶跳過繁瑣的軟件安裝及環(huán)境配置，簡化分析流程，提高了數(shù)據(jù)處理的效率和處理的可行性[8]。傳統(tǒng)的生物信息分析模式是單任務(wù)處理模式，對于大量測序數(shù)據(jù)來說，碎片化分析工作過于繁瑣，構(gòu)建本系統(tǒng)來分析則大大減輕了復(fù)雜程度。在大數(shù)據(jù)時代，高效分析顯得尤為重要。本系統(tǒng)處理任務(wù)時，只需要將完整的單細(xì)胞轉(zhuǎn)錄組測序數(shù)據(jù)上傳至指定工作目錄，運(yùn)行Python 自動化腳本即可依照預(yù)設(shè)流程，完成從數(shù)據(jù)讀取質(zhì)控到分析計(jì)算結(jié)束的工作，并通過SFTP 文件服務(wù)器下載所需的結(jié)果文件。應(yīng)用本系統(tǒng)可以降低分析人員使用Linux 環(huán)境的難度，也避免了多用戶在服務(wù)器直接操作帶來的宕機(jī)風(fēng)險(xiǎn)和不安全因素，可以最大化服務(wù)器的工作效率和穩(wěn)定性[9]。

本研究開發(fā)的單細(xì)胞轉(zhuǎn)錄組病毒信息分析系統(tǒng)仍存在一定的局限性。從理論方法層面上看，測試所使用的數(shù)據(jù)集數(shù)量比較少，還有待增加。系統(tǒng)最大負(fù)荷運(yùn)算能力尚未評估，系統(tǒng)滿負(fù)荷穩(wěn)定性尚未評估，此方面后續(xù)研究需改進(jìn)[10]；使用層面上，需要使用者具備一定的Linux 使用基礎(chǔ)以及生物信息學(xué)知識，個別參數(shù)針對不同類型來源數(shù)據(jù)需調(diào)整，還未開發(fā)人性化和諧美觀的使用交互界面，系統(tǒng)美觀性還存在不足。