少量細胞輸入的自制轉錄組測序文庫構建試劑評測

丁 蕾 ，高彩霞 ，劉兆遠 ，陳 磊

1.上海交通大學基礎醫學院免疫學與微生物學系，上海 200025；2.上海交通大學醫學院，上海市免疫學研究所，上海 200025

伴隨著高通量測序技術的迅猛發展和測序成本的降低，轉錄組測序借助其靈敏度高、檢測范圍廣、提供轉錄組信息全面等獨特的優勢在轉錄組學研究領域逐漸占據主導地位[1]，是目前深入研究轉錄組復雜性的強大工具。

針對轉錄組測序的生物信息學分析主要包括以下幾個步驟：數據質量控制、序列比對和基因表達分析。其中數據質量控制主要包括評估建庫試劑在實驗中捕獲轉錄本的能力和所測序列在基因覆蓋度的均勻性。基因表達分析主要包括差異表達基因的鑒定和差異表達基因富集分析，如KEGG（Kyoto encyclopedia of genes and genomes）通路富集分析，這是轉錄組測序最重要的分析結果。基因表達分析結果可用于更深入的研究，包括轉錄本結構研究（如基因的可變剪切）、轉錄本變異研究（如基因融合、單核苷酸突變）、差異基因表達水平的比較，甚至包括全新轉錄本或稀有轉錄本的發現[2]。尤其對于如腫瘤干細胞等來源極為有限的生物樣品，轉錄組測序技術的優勢明顯[3]，要求的起始樣品量要比芯片技術少得多且技術重復性好。SMART（switching mechanism at 5′ end of RNA template）技術經優化后，選擇使用鎖核酸（locked nucleic acid，LNA）、更高濃度的氯化鎂以及甜菜堿的新方案，使得起始僅使用1～1 000個細胞或10 pg～10 ng的總RNA就可獲得高質量的測序文庫，而且還使得序列在基因的覆蓋度上獲得很大改善[4]，能夠更好地實現單核苷酸變異等的檢測。基于該技術的商業化試劑盒有多種，均價格昂貴，其中普遍使用的是TaKaRa Bio公司名為SMART-Seq v4 Ultra Low Input RNA Kit的試劑盒（下文簡稱為TaKaRa試劑）。本研究利用SMART技術自制相對價格更低的建庫試劑（下文簡稱為DIY試劑）用于少量細胞輸入的轉錄組測序文庫構建實驗，并通過生物信息學分析，從多個方面比較DIY試劑和TaKaRa試劑對轉錄組測序結果的影響，從而驗證用DIY試劑替代TaKaRa試劑的可行性，以降低轉錄組測序文庫構建的成本。

1 材料與方法

1.1 材料

1.1.1 實驗細胞和動物 小鼠腹腔巨噬細胞，分別來自于4只8周齡雌性SPF級C57BL/6小鼠。小鼠飼養于上海交通大學醫學院實驗動物科學部屏障環境內。生產許可證號為SCXK（滬）2018-0007，使用許可證號為SYXK（滬）2018-0027。

1.1.2 主要試劑 TaKaRa試劑：SMART-Seq v4 Ultra Low Input RNA Kit，購于TaKaRa Bio公司。DIY試劑：不含鈣、鎂的磷酸緩沖液DPBS（Dulbecco′s phosphate buffered saline），1×TrypLE Express 酶（無酚紅），購于Gibco公司；RNaseZap，購于Ambion公司；DNA去除試劑（DNAOFF），購于TaKaRa Bio公司；2′-脫氧核苷酸-5′-三磷酸混合物（10 mmol/L dNTP mix），購 于Fermentas公 司；5×First-strand Buffer（250 mmol/L Tris-HCl，pH 8.3，室溫；375 mmol/L 氯化鉀；15 mmol/L 氯化鎂；二硫蘇糖醇）、Superscript Ⅱ Reverse Transcriptase，購于Invitrogen公司；重組RNase 抑制劑，購于Clonetech公司；甜菜堿、氯化鎂（無水），購于Sigma-Aldrich公司。其他試劑：蒸餾水，購于Gibco公司；2×KAPA HiFi HotStart ReadyMix，購于KAPA Biosystems公司；Agencourt Ampure XP beads，購于Beckman Coulter公司；溴化乙錠（ethidium bromide，EB，10 mmol/L Tris-HCl，pH 8.5），購于Qiagen公司；雙端測序TruSeq Dual-index Sequencing Primer試劑盒、Nextera XT DNA Sample Preparation試劑盒、Nextera XT 24-index試劑盒、Adapter oligos，購于Illumina公司；99.5%乙醇，購于Sigma-Aldrich公司；TRIzol 試劑、Glycoblue，購于Invitrogen公司；氯仿、異丙醇，購于國藥試劑公司；4%巰基乙酸鹽肉湯，購于BD公司。

1.1.3 主要儀器 Agilent 2100 生物分析儀系統，購于安捷倫科技有限公司；Illumina NextSeq 500儀器，購于Illumina公司；移液槍、低溫超速離心機、高速臺式離心機，購于Eppendorf公司；EasyCycler 96 PCR儀，購于Analytik Jena AG公司。

1.1.4 主要軟件 使用基于JAVA 8 的FastQC 0.11.9對高通量測序數據進行質量評估；使用基于Python 3.7 的RSeQC 3.0.1 評估高通量測序尤其是轉錄組測序數據的覆蓋均勻性；使用基于C++ 的fastp 0.20.0對原始測序數據進行隨機抽樣等操作；使用HISAT2 （hierarchical indexing for spliced alignment of transcripts 2） 2.1.0進行有參考基因組的序列比對；使用基于C 語言的SAMtools 1.4處理經HISAT2 等軟件比對后生成的sam/bam 文件；基于Python 3.7 的HTSeq 0.11.3處理高通量測序數據；使用基于R 3.6.1的DESeq2、edgeR、ggplot2軟件包進行基因差異分析以及數據可視化分析等。

1.2 實驗方法

1.2.1 腹腔巨噬細胞的誘導 將4只小鼠隨機均分為2組，再分別編號為M1、M2、M3、M4，其中M1、M2號小鼠作為對照組不進行任何處理，M3、M4號小鼠作為實驗組向其腹腔內注射1 mL濃度為4%的巰基乙酸鹽肉湯以誘導巨噬細胞[5]。經72 h后，利用流式細胞分選技術分離腹腔巨噬細胞。

1.2.2 腹腔巨噬細胞RNA的提取 腹腔巨噬細胞離心5 min，4 ℃，500×g；向每個樣品中加入500 μL TRIzol，混合均勻，室溫孵育5 min；置于-80 ℃至少1 h；向每個樣品中添加適量的氯仿（TRIzol與氯仿的體積比為5:1），劇烈搖動約15 s，室溫孵育2～3 min；離心15 min，4 ℃，12 000×g；將水相轉移至新管中（TRIzol與氯仿的體積比為2:1），此后在冰上工作；向每個樣品中添加1 μL共沉淀劑GlycoBlue；在每個樣品中加入適量的異丙醇（TRIzol與異丙醇的體積比為2:1），混合均勻，-80 ℃孵育1 h以上；離心20 min，4 ℃，12 000×g；用適量的75%預冷過的乙醇洗滌（TRIzol與異丙醇的體積比為1:1）；離心15 min，4 ℃，7 400×g，風干以盡可能地除去乙醇；在12.5～15.0 μL無核酸酶的水中將每個RNA沉淀重懸；室溫孵育2～3 min，充分混合并置于冰上；使用Qubit儀器定量RNA；保存RNA于-80℃條件下。

1.2.3 cDNA的構建及二代測序 將編號為M2和M3的小鼠的腹腔巨噬細胞提取的RNA各均分成2份，用DIY試劑和TaKaRa試劑分別進行cDNA的構建，所得4個樣品的編號為C_M2_DIY、C_M2_TaKaRa、T_M3_DIY和T_M3_TaKaRa（C，對照組；T，實驗組）。另外，從編號為M1和M4的小鼠的腹腔巨噬細胞提取的RNA皆用DIY試劑進行cDNA的構建，所得2個樣品的編號為C_M1_DIY、T_M4_DIY。使用Agilent 2100 生物分析儀系統對6個樣品的cDNA進行質檢。在完成轉錄組測序文庫構建后，對上述6個樣品利用Illumina NextSeq 500平臺進行雙端測序，測序模式為2×75 bp。

1.2.4 測序數據質量控制和序列比對 使用FastQC軟件對測序數據進行初步質量控制，輸出結果會給每個堿基一個相應的質量評分，用于衡量測序精確度。一個給定堿基的測序質量評分Q定義為：Q=-10×lge，其中e為預計堿基檢出不正確的概率。故堿基的質量評分與堿基檢出精確度存在一定關系，即Q較高表示出錯的概率較小。其中，質量控制結果中的重要指標Q30（%）代表堿基質量評分≥30的堿基的數量占全部堿基數的百分比。一般而言，Q30＞85%時，測序數據質量合格。使用fastp軟件對原始的序列文件進行隨機抽樣，使用HISAT2軟件[6]進行有參考基因組的序列比對，并使用SAMtools軟件[7]處理比對結果，再用HTSeq軟件[8]統計比對到參考基因組上基因區間內的讀序數目，計算不同測序數據量下可檢測的轉錄本數，比較DIY試劑和TaKaRa試劑捕獲轉錄本的能力。使用RSeQC軟件[9]對測序數據進行質量控制，獲得測序數據在基因上的覆蓋均勻性。從數據質量以及序列比對的層面探究2種試劑對轉錄組測序結果的影響。

1.2.5 基因表達分析 經過HTSeq軟件完成基因表達量的計算后，使用R語言軟件包DESeq2[10]和edgeR[11]進行樣品的相似性分析以及實驗組與對照組之間的差異表達基因分析[12]，選擇符合|log2fold change|＞1、P＜0.005的基因為差異表達基因，探究2種試劑處理樣品的差異表達基因鑒定結果是否一致。log2fold change表示以2為底、基因表達量差異的對數，代表實驗組中的基因表達量與對照組中的基因表達量的差異倍數。若log2fold change的數值為正，則反映該基因在實驗組中的表達量比在對照組中的表達量高；相反若數值為負，則反映該基因在實驗組中的表達量比在對照組中的表達量低。

1.2.6 差異表達基因通路富集分析 使用在線富集分析工具DAVID（database for annotation, visualization, and integrated discovery）6.8[13]對差異表達基因進行通路富集分析，并使用R語言軟件包ggplot2[14]進行可視化分析。通路富集分析方法以通路為單位，以物種已知的全部基因為背景，通過Fisher精確檢驗來分析計算各個通路基因富集度的顯著性水平，從而確定受到顯著影響的代謝和信號轉導途徑，探究2種試劑處理情況下差異表達基因通路富集分析的結果是否一致。以錯誤發現率（false discovery rate，FDR）表示富集的顯著性水平。FDR表示以10為底、P值的對數，FDR數值越大，代表在統計學上該通路的影響或變化越顯著。

2 結果

2.1 2種試劑處理樣品的數據質量比較

6個樣品的cDNA片段均分布在200～3 000 bp范圍內，主峰均在1 000 bp左右，cDNA質量良好。使用FastQC軟件對測序原始數據進行質量控制，結果如表1所示。Q30均＞91%，表明數據質量良好。經HISAT2軟件比對后，發現比對到基因組的序列占比均＞90%，捕獲的轉錄本數也基本相近。

表1 樣品基本信息及數據質量控制Tab 1 Sample basic information and data quality control

2.2 2種試劑捕獲轉錄本的能力比較

在一定范圍內，轉錄組測序結果顯示捕獲到的轉錄本數會隨著測序數據量的增加而增加，并最終趨于飽和。為了比較DIY試劑和TaKaRa試劑捕獲轉錄本的能力，使用fastp軟件對原始的序列文件進行隨機抽樣，獲得具有0.5、1、1.5、2、2.5、3、4、5、6、7 M的序列子集文件，再用HTSeq軟件計算不同測序數據量下可檢測的轉錄本數。如圖1所示：可檢測的轉錄本數隨著測序數據量的增加而增加；測序數據量一定時，對照組（C_M1_DIY、C_M2_DIY、C_M2_TaKaRa）和 實 驗 組（T_M3_DIY、T_M4_DIY、T_M3_TaKaRa）可檢測的轉錄本數差異明顯，而不同試劑處理的樣品各自在對照組（C_M2_DIYvsC_M2_TaKaRa）和實驗組（T_M3_DIYvsT_M3_TaKaRa）內無明顯差異。

圖1 2種試劑捕獲轉錄本的能力比較Fig 1 Comparison of the capability in transcripts capture of DIY reagent and TaKaRa reagent

2.3 2種試劑處理的樣品基因覆蓋度比較

使用RSeQC軟件計算所測序列在基因上的覆蓋度，結果如圖2所示。DIY試劑和TaKaRa試劑處理的樣品基因覆蓋度都較均勻、無偏差，且兩者的一致性較高。

圖2 2種試劑處理的樣品基因覆蓋度Fig 2 Gene coverage plot of the samples treated with the two reagents

2.4 主成分分析

分別使用R軟件的DESeq2包和ggplot2包對樣品進行主成分分析（principal component analysis，PCA）以及可視化分析。圖3清晰地體現了樣品的聚類情況，橫坐標PC1反映實驗處理導致的樣品間差異，即巰基乙酸鹽肉湯對腹腔巨噬細胞持續72 h的誘導作用，存在91%的差異性；縱坐標PC2反映2種試劑導致的樣本間差異，該差異性僅為8%，說明DIY試劑和TaKaRa試劑對樣品的影響遠小于實驗處理導致的基因表達譜差異。而且，DIY試劑處理的對照組樣品（C_M1_DIYvsC_M2_DIY）和實驗組樣品（T_M3_DIYvsT_M4_DIY）的組內差異均較小，表示DIY試劑的實驗重復性好且結果穩定。

圖3 PCA體現的樣品聚類情況Fig 3 Sample clustering reflected by PCA

2.5 2種試劑處理樣品的差異表達基因鑒定結果

將6個樣品分為2組：第1組為DIY試劑處理的樣品，包括C_M1_DIY、C_M2_DIY、T_M3_DIY、T_M4_DIY；第2組為TaKaRa試劑處理的樣品，包括C_M2_TaKaRa和T_M3_TaKaRa。使用edgeR包對2組數據就對照組和實驗組之間的差異表達基因進行鑒定，最終DIY試劑處理組樣品共有3 877個差異表達的基因，TaKaRa試劑處理組樣品共有3 855個差異表達的基因，重合的基因數為2 737個。結果如圖4所示，從2組結果中各取顯著性差異最大的前1 000個基因（按照P值排序），重合的基因數有748個，重合率達75%，并且相同基因在2組數據中的|log2fold change|值比較接近。

圖4 2種試劑處理樣品的差異表達基因鑒定結果Fig 4 Differential gene expression analysis results of the samples treated with the two reagents

2.6 2種試劑處理樣品的基因通路富集分析

根據上述按照所用試劑將樣品分為2組的策略，各自選取差異表達分析結果中最顯著的1 000個基因（按照P值排序），使用DAVID在線分析工具進行KEGG通路富集分析，結果顯示2組差異基因中富集程度顯著的通路一致性較高，重合率達80%。使用R語言的ggplot2軟件包對各自重合的前10個通路進行可視化，結果如圖5所示。橫坐標數值為正數代表在實驗組中上調基因富集出來的結果，為負值則代表在實驗組中下調基因富集出來的結果；數值的大小則體現通路富集的顯著性水平。結果表明，DIY試劑處理組和TaKaRa試劑處理組的基因富集通路的顯著差異情況較相近。

圖5 2種試劑處理樣品的基因通路富集分析結果Fig 5 Results of KEGG pathway analysis of the samples treated with the two reagents

3 討論

轉錄組測序的研究對象為特定細胞在某一功能狀態下能轉錄出來的所有RNA的總和，可以用于基因功能以及基因結構的研究[15]。轉錄組測序結果可以揭示特定生物學過程以及疾病發生過程中的分子機制，因此該技術被廣泛應用于生物學、基礎醫學、臨床診斷和藥物研發等多個領域[16]，并成為一個識別細胞類型、鑒定標記基因、識別信號通路和研究調控機制的重要工具。轉錄組測序技術可為不同個體和疾病的組織帶來更好的生物靶標識別和藥物標靶[17]，從而有助于進一步發展精準醫學。

本研究使用DIY試劑和TaKaRa試劑分別進行轉錄組測序文庫構建，并通過生物信息學分析，從不同的方面來驗證DIY試劑代替TaKaRa試劑的可行性。結果發現，2種試劑處理樣品得到的數據質量皆良好，比對到基因組的序列占比也基本相等；在相同的測序深度下，2種試劑捕獲的轉錄本數較相近；2種試劑處理的樣品所測序列在基因上的覆蓋度均勻，且一致性較高。此外，分析結果顯示DIY試劑在僅有少量細胞輸入的情況下，能夠捕獲相對較多的轉錄本數量，表明該試劑使得來源極為有限的生物樣品分析成為可能。所以從數據質量、捕獲轉錄本能力、基因覆蓋度均勻性等方面來看，DIY試劑可以很好地替代TaKaRa試劑進行少量細胞輸入的轉錄組測序文庫構建。差異表達基因的鑒定和差異基因通路富集分析是轉錄組測序較重要的分析結果[18]，所以有必要比較這2個分析結果在2種不同試劑處理的情況下是否一致。結果顯示，2組分析結果相近，無明顯差別。從構建文庫的成本方面來看，DIY試劑的優勢更加明顯：若用于建庫的起始細胞量為100個細胞，TaKaRa試劑建庫實驗所花費用高達6萬元，DIY試劑則僅需約5 700元。然而在進行一項科學研究時，動輒需要對上千萬乃至更大數量級的細胞進行文庫構建及測序，測序的費用更是昂貴，所以使用DIY試劑代替此類商業化試劑不僅能在保證數據質量良好、差異基因分析等結果可靠的前提下降低約90%的建庫成本，還能根據具體實驗進行酶反應條件的調整以及試劑使用的優化等。因此綜合多方面的評測，DIY試劑可以替代昂貴的商業化試劑用于少量細胞輸入的轉錄組測序文庫構建。此外，SMART技術不僅可用于少量細胞輸入的轉錄組測序文庫構建，也可對單個細胞進行建庫并測序，用以在免疫、神經等復雜系統中研究細胞高異質性這一類問題；該技術不僅為測序技術提供了新的發展方向，還有望在單細胞水平上進行遺傳變異的檢測、腫瘤的診斷及免疫治療等研究[19]。因此，如何優化DIY試劑并將其用于單細胞的轉錄組測序文庫構建還亟待解決。

綜上所述，本研究比較了DIY試劑與TaKaRa試劑對少量細胞的轉錄組測序文庫構建以及轉錄組測序結果，表明DIY試劑可以替代昂貴的商業化試劑進行轉錄組測序文庫構建。該研究將使得基于SMART技術的轉錄組測序得到真正的廣泛應用，將為生物學以及基礎醫學等各個研究領域提供良好的技術支持。