轉(zhuǎn)錄組測序的研究對象為特定細(xì)胞在某一功能狀態(tài)下所能轉(zhuǎn)錄出來的所有 RNA 的總和，包括 mRNA 和非編碼 RNA，相對于傳統(tǒng)的芯片雜交平臺，轉(zhuǎn)錄組測序無需預(yù)先針對已知序列設(shè)計探針，即可對任意物種的整體轉(zhuǎn)錄活動進(jìn)行檢測，提供更準(zhǔn)確的數(shù)字化信號，更高的檢測通量以及更廣泛的檢測范圍，是目前深入研究轉(zhuǎn)錄組復(fù)雜性的強(qiáng)大工具，基于高通量測序平臺的轉(zhuǎn)錄組測序技術(shù)能夠全面獲得物種特定組織或器官的轉(zhuǎn)錄本信息，從而進(jìn)行基因表達(dá)水平研究、新轉(zhuǎn)錄本發(fā)現(xiàn)研究、轉(zhuǎn)錄本結(jié)構(gòu)變異研究等

• RNA-seq 概述

• 版本信息

• 一些名詞解釋

• 大致流程

• 軟件安裝

• 數(shù)據(jù)獲取與預(yù)處理

• 測序數(shù)據(jù)下載

• 參考基因組及注釋文件

• FastQc 質(zhì)控

• 常用參數(shù)

• 進(jìn)行質(zhì)量檢測

• Trimmomatic 過濾低質(zhì)量序列

• 常用參數(shù)

• 過濾低質(zhì)量序列

• hisat2 比對

• 建立索引

• 常用參數(shù)

• 進(jìn)行比對

• samtools 排序壓縮

• 常用命令

• 排序壓縮

• featureCounts 生成基因計數(shù)表

• 常用參數(shù)

• 計數(shù)統(tǒng)計

• 參考

RNA-seq 概述

RNA-seq 是研究轉(zhuǎn)錄組應(yīng)用最廣泛，也是最重要的技術(shù)之一，RNA-seq 分析內(nèi)容包括序列比對、轉(zhuǎn)錄本拼接、表達(dá)定量、差異分析、融合基因檢測、可變剪接、RNA 編輯和突變檢測等，具體流程和常用工具如下圖所示，通常的分析不一定需要走完全部流程，按需進(jìn)行，某些步驟可以跳過、簡化等

RNA-seq 中最常用的分析方法就是找出差異表達(dá)基因 (Differential gene expression, DEG)，在實驗室中，標(biāo)準(zhǔn)流程就分為三步：

• step1: 構(gòu)建測序文庫，包括提取 RNA, 富集 mRNA 或清除核糖體 RNA, 合成 cDNA, 加上接頭

• step2: 在高通量測序平臺（通常為 Illumina） 上對文庫進(jìn)行測序，每個樣本的測序深度為 10-30M 讀長

• step3: 數(shù)據(jù)分析，具體而言：對測序得到的讀長進(jìn)行比對或組裝到轉(zhuǎn)錄本上；對覆蓋到每個基因區(qū)域的讀長進(jìn)行計數(shù)；根據(jù)統(tǒng)計模型鑒定不同樣本間差異表達(dá)的基因，（這種分析過程是比較傳統(tǒng)的方法）

版本信息

• 成都理工大學(xué)超算平臺 Red Hat 4.8.5-36

• conda 4.10.3

• 測序數(shù)據(jù)來自 NCBI ACC=SRR25909836

• 參考基因組和注釋文件來自 Ensembl Plants

• 物種：Oryza sativa Japonica

一些名詞解釋

• adapter：接頭，為一段已知的短核苷酸序列，用于鏈接未知的目標(biāo)測序片段

• index：幾個堿基組成的寡核苷酸鏈，用于在混合測序時，區(qū)分不同樣本

• insert：待測序的目標(biāo)序列，位于兩個 adapter 之間

• reads：在測序過程中，從樣本中識別出來的 DNA 或 RNA 序列

• Transcriptome: 包含所有 RNA 分子的集合，可以用來量化每個基因的表達(dá)水平

• Expression levels 或 Coverage: 表示特定基因的 RNA 數(shù)量，通常用來衡量基因的表達(dá)水平

• Differential expression analysis：用來比較不同樣本或不同處理組之間的基因表達(dá)差異

• Gene annotation：一個包含有關(guān)每個基因的基本信息的數(shù)據(jù)庫，如其位置，功能等

• Splice variants：一個基因能夠通過不同的剪接方式生成不同的 RNA，這些不同的 RNA 就叫做剪接異構(gòu)體

• FPKM (Fragments Per Kilobase of transcript per Million mapped reads) 和 TPM (Transcripts Per Kilobase Million)：用作基因表達(dá)水平的衡量標(biāo)準(zhǔn)

• Quality control：確保數(shù)據(jù)質(zhì)量，清理模糊不清或者質(zhì)量低下的 reads

• Genome：一個生物體的所有遺傳信息

大致流程

軟件安裝

只需要使用 conda 就可以安裝所有需要的軟件，主要使用的軟件有以下一些

• sra-tools：快速下載 NCBI SRA 數(shù)據(jù)

• FastQc：測序數(shù)據(jù)質(zhì)量檢測與控制

• Trimmomatic：過濾低質(zhì)量序列

• hisat2：轉(zhuǎn)錄組數(shù)據(jù)的比對

• samtools：對 hisat2 比對的結(jié)果進(jìn)行排序和壓縮

• featureCounts：對基因的信息進(jìn)行計數(shù)統(tǒng)計

先創(chuàng)建一個虛擬環(huán)境，要注意的是最新版的 samtools 需要 3.9 及以上環(huán)境

conda create -n rna-seq python=3.10

安裝所有需要的軟件

conda install fastqc trimmomatic hisat2 subread #subread 即 featureCounts

SRA Toolkit 和 samtools 安裝時有坑，先搜索 sra-tools 和 samtools 的版本，然后安裝指定版本

conda search sra-tools
conda install sra-toolst=3.0.7

conda search samtools
conda install samtools=1.9

數(shù)據(jù)獲取與預(yù)處理

測序數(shù)據(jù)下載

先在 NCBI 的 SRA 數(shù)據(jù)庫搜索感興趣的物種

選擇符合自己要求的文章，找到下面 Runs 這里，點擊 SRR 開頭的編號

查看數(shù)據(jù)是否符合要求

文件是***_1.fq.gz、***_2.fq.gz 這種是雙端測序數(shù)據(jù)，我們需要這種雙端測序的數(shù)據(jù)來進(jìn)行 RNA-seq 分析

如果數(shù)據(jù)是雙端測序的，那么就復(fù)制 SRR 編號，使用 sratools 下載，例如：

prefetch SRR8956151

批量下載需要先建立一個 txt 文件，將 SRR 編號寫進(jìn)去，例如：

SRR5830630
SRR5830631
SRR5830632
SRR5830633
SRR5830634

然后使用下面的命令下載

prefetch --option-file SRR_Acc_List.txt

由于數(shù)據(jù)比較大，可以使用 nohup 命令掛在后臺下載

nohup prefetch --option-file SRR_Acc_List.txt

剛剛下載好的數(shù)據(jù)是 sra 格式的，使用 sratools 將其拆分

fastq-dump --gzip --split-3 SRR25909836.sra

• –gzip 是將拆分的 fastq 文件壓縮歸檔為 gz 格式

• –split-3 是將文件拆分為正向序列和逆向序列

如果數(shù)據(jù)比較多，就寫一個 bash 腳本

#!/bin/bash
mkdir SRR
mv ./SRR*/*.sra ./SRR
cd SRR
nohup fastq-dump --gzip --split-3 SRR*.sra

參考基因組及注釋文件

植物的我一般在 Ensembl Plants 下載，用 wget 或 curl 都可以，內(nèi)存不大

wget https://ftp.ensemblgenomes./pub/plants/release-57/fasta/oryza_sativa/dna/Oryza_sativa.IRGSP-1.0.dna.toplevel.fa.gz

wget https://ftp.ensemblgenomes./pub/plants/release-57/gff3/oryza_sativa/Oryza_sativa.IRGSP-1.0.57.gff3.gz

然后解壓

gzip -d Oryza_sativa.IRGSP-1.0.57.gff3.gz
gzip -d Oryza_sativa.IRGSP-1.0.dna.toplevel.fa.gz

為了方便，我把兩個文件分別重命名為 oryza_sativa.fa 和 oryza_sativa.gff3

mv Oryza_sativa.IRGSP-1.0.dna.toplevel.fa oryza_sativa.fa
mv Oryza_sativa.IRGSP-1.0.57.gff3 oryza_sativa.gff3

FastQc 質(zhì)控

常用參數(shù)

fastqc [-o output dir] [–(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN

• -o –outdir：FastQC 生成的報告文件的儲存路徑，生成的報告的文件名是根據(jù)輸入來定的

• –extract：生成的報告默認(rèn)會打包成 1 個壓縮文件，使用這個參數(shù)是讓程序不打包

• -t –threads：選擇程序運(yùn)行的線程數(shù)，每個線程會占用 250MB 內(nèi)存，越多越快咯

• -c –contaminants：污染物選項，輸入的是一個文件，格式是 Name [Tab] Sequence，里面是可能的污染序列，如果有這個選項，F(xiàn)astQC 會在計算時候評估污染的情況，并在統(tǒng)計的時候進(jìn)行分析，一般用不到

• -a –adapters：也是輸入一個文件，文件的格式 Name [Tab] Sequence，儲存的是測序的 adpater 序列信息，如果不輸入，目前版本的 FastQC 就按照通用引物來評估序列時候有 adapter 的殘留

• -q –quiet：安靜運(yùn)行模式，一般不選這個選項的時候，程序會實時報告運(yùn)行的狀況

進(jìn)行質(zhì)量檢測

這里我們使用 fastqc [文件名] 即可

fastqc SRR25909836_1.fastq.gz

當(dāng)然，數(shù)據(jù)比較多的時候還是掛在后臺批處理然后等著就行

nohup fastqc SRR*.fastq.gz

程序運(yùn)行完成后會輸出一堆 html 文件和 zip 壓縮包，html 是網(wǎng)頁版報告，zip 是本地寶報告，下載到本地用瀏覽器打開就可以看到質(zhì)量檢測報告了

左側(cè) Summary 部分就是整個報告的目錄，整個報告分成若干個部分

• 合格：綠色的√

• 警告：黃色的！

• 不合格：紅色的×

我們一般比較關(guān)心的是下面幾個部分

• Basic Statistics: 對數(shù)據(jù)量的概覽

• Per base sequence quality：reads 每個位置測序質(zhì)量最直接的展示

• Per sequence quality scores：總體 reads 測序質(zhì)量趨勢

• Per base sequence content：ATGC 含量估計測序是否存在偏差

• Sequence Duplication Levels：影響測序的因素太多，查看是否存在污染，數(shù)據(jù)處理時是否需要去冗余；現(xiàn)在數(shù)據(jù)量都可以滿足需求，因此前期數(shù)據(jù)處理時，盡量高標(biāo)準(zhǔn)，嚴(yán)格質(zhì)控

Basic Statistics

每個位置的堿基的測序質(zhì)量

• Encoding：指測序平臺的版本和相應(yīng)的編碼版本號

• Total Sequences：總的 reads 數(shù)

• Sequence length：測序的長度

• %GC：是我們需要重點關(guān)注的一個指標(biāo)，這個值表示的是整體序列中的 GC 含量，這個數(shù)值一般是物種特異的，比如人類細(xì)胞就是 42%左右，如果測序原始數(shù)據(jù)的 GC 含量遠(yuǎn)遠(yuǎn)偏離這個比例，說明測序數(shù)據(jù)存在一定偏好性，如果直接用測序數(shù)據(jù)，會影響后續(xù)的 CNV 和變異檢測的分析

Per base sequence quality

• 橫軸代表位置（第 1 到 150 個堿基）

• 縱軸代表 quality

• 紅色表示中位數(shù)

• 黃色是 25%-75%區(qū)間

• 觸須是 10%-90%區(qū)間

• 藍(lán)線是平均數(shù)

• Warning，如果任何堿基質(zhì)量低于 10, 或者是任何中位數(shù)低于 25

• Failure，如果任何堿基質(zhì)量低于 5, 或者是任何中位數(shù)低于 20

一般要求此圖中，所有位置的 10%分位數(shù)大于 20，否則切除 20 以下的堿基，從而保證后續(xù)分析的正確性

Per tile sequence quality

每個 tile 測序的測序質(zhì)量

• tile：每一次測序熒光掃描的最小單位

• 橫軸代表 101 個堿基的位置

• 縱軸是 tail 的 Index 編號

檢查 reads 中每一個堿基位置在不同的測序小孔之間的偏離度，藍(lán)色表示低于平均偏離度，偏離度小，質(zhì)量好；越紅表示偏離平均質(zhì)量越多，質(zhì)量也越差，如果出現(xiàn)質(zhì)量問題可能是短暫的，如有氣泡產(chǎn)生，也可能是長期的，如在某一小孔中存在殘骸，問題不大，可以看到我的這個數(shù)據(jù)幾乎沒有瑕疵

Per sequence quality scores

每條序列質(zhì)量得分的分布情況

• 橫軸是平均質(zhì)量得分

• 縱軸收到該平均質(zhì)量得分的 reads 數(shù)目

• 當(dāng)測序質(zhì)量峰值小于 27（錯誤率 0.2%）時報"WARN"

• 當(dāng)峰值小于 20（錯誤率 1%）時報"FAIL"

假如我測的 1 條序列長度為 101bp，那么這 101 個位置每個位置 Q 值的平均值就是這條 reads 的質(zhì)量值，我這個數(shù)據(jù)的質(zhì)量不錯，reads 大都集中在高分上，縱軸數(shù)值越大，該序列測序錯誤的可能就越小

Per base sequence content

統(tǒng)計 reads 每個位置 ATCG 四種堿基的分布

• 橫軸是 1 - 101 bp；縱軸是百分比

• 圖中四條線代表 A T C G 在每個位置平均含量

• 理論上來說，A 和 T 應(yīng)該相等，G 和 C 應(yīng)該相等，但是一般測序的時候，剛開始測序儀狀態(tài)不穩(wěn)定，很可能出現(xiàn)上圖的情況，像這種情況，即使測序的得分很高，也需要 cut 開始部分的序列信息

Per sequence GC content

序列平均 GC 含量分布

• 橫軸是百分比

• 縱軸是每條序列 GC 含量對應(yīng)的數(shù)量

• 藍(lán)色的線是程序根據(jù)經(jīng)驗分布給出的理論值

• 紅色是真實值

藍(lán)紅色線應(yīng)該比較接近才比較好，當(dāng)紅色的線出現(xiàn)雙峰，很有可能是混入了其他物種的 DNA 序列，比如我這張圖

Per base N content

每個位置無法檢測的值的比例，當(dāng)測序儀無法確定是何種堿基時，用 N 表示

• 當(dāng)任意位置的 N 的比例超過 5%，報"WARN"

• 當(dāng)任意位置的 N 的比例超過 20%，報"FAIL"

正常情況下，N 的比例是很小的，所以圖上常?？吹揭粭l直線，但放大 Y 軸之后會發(fā)現(xiàn)還是有 N 的存在，這不算問題，當(dāng) Y 軸在 0%-100%的范圍內(nèi)也能看到“鼓包”時，說明測序系統(tǒng)出了問題

Sequence Length Distribution

序列測序長度分布

每次測序儀測出來的長度在理論上應(yīng)該是完全相等的，但是總會有一些偏差，當(dāng)測序的長度不同時，如果很嚴(yán)重，則表明測序儀在此次測序過程中產(chǎn)生的數(shù)據(jù)不可信，比如我的這個圖中，150pb 是最主要的，其他的幾乎沒有，所以數(shù)據(jù)的質(zhì)量還是比較高的

Sequence Duplication Levels

統(tǒng)計 reads 重復(fù)水平

• 橫坐標(biāo)是重復(fù)的次數(shù)

• 縱坐標(biāo)是 duplicated reads 占 unique reads 總數(shù)百分比

• 當(dāng)非 unique 的 reads 占總數(shù)的比例大于 20%時，報"WARN"

• 當(dāng)非 unique 的 reads 占總數(shù)的比例大于 50%時，報"FAIL"

測序本身就會產(chǎn)生重復(fù) reads, 測序深度越高，reads 重復(fù)數(shù)越大；如果重復(fù)出現(xiàn)峰值，就提示可能存在偏差（如建庫過程中的 PCR duplication）

fastqc 抽取 reads 文件前 200,000 條 reads 統(tǒng)計其重復(fù)情況，重復(fù)數(shù)目大于等于 10 的 reads 被合并統(tǒng)計，這也是為什么我們看到上圖的中間那里略有上揚(yáng)，大于 75bp 的 reads 只取 50bp 進(jìn)行比較，由于 reads 越長錯誤率越高，所以其重復(fù)程度仍有可能被低估

Overrepresented sequences

過度重復(fù)出現(xiàn)的序列的統(tǒng)計信息，上圖中沒有

Adapter Content

衡量的是序列中兩端 adapter 的情況

如果在當(dāng)時 fastqc 分析的時候-a 選項沒有內(nèi)容，則默認(rèn)使用圖例中的四種通用 adapter 序列進(jìn)行統(tǒng)計

上圖中 adapter 都已經(jīng)去除，如果有 adapter 序列沒有去除干凈的情況，在后續(xù)分析的時候需要先使用 cutadapt 等軟件進(jìn)行去接頭

Trimmomatic 過濾低質(zhì)量序列

常用參數(shù)

• PE：雙端測序文件過濾

• SE：單端測序文件過濾

• -baseout：輸出目錄

• -threads：線程數(shù)，最大是 CPU 核數(shù)

• -trimlog：生成日志名

• -quiet：靜默模式

• ILLUMINACLIP：從 reads 中剪切 adapter 和其他 Illumina 特定序列

• SLIDINGWINDOW：執(zhí)行滑動窗口修剪，一旦窗口內(nèi)的平均質(zhì)量低于閾值，則切割

• LEADING：如果低于閾值質(zhì)量，則在 reads 起始處剪切堿基

• TRAILING：如果低于閾值質(zhì)量，則在 reads 末尾處剪切堿基

• CROP：將 reads 從末尾切割為指定長度

• HEADCROP：從 reads 剪切后低于指定長度，則刪除

• MINLEN：如果 reads 低于指定長度，則刪除

• TOPHRED33：將質(zhì)量得分轉(zhuǎn)換為 Phred-33

• TOPHRED64：將質(zhì)量得分轉(zhuǎn)換為 Phred-64

過濾低質(zhì)量序列

我使用的是下面的命令，需要根據(jù)自己的文件進(jìn)行調(diào)整

trimmomatic PE -threads 1 -phred33 SRR25909836_1.fastq.gz SRR25909836_2.fastq.gz -summary oryza_sativa.summary -baseout SRR25909836.fastq.gz LEADING:3 TRAILING:3 SLIDINGWINDOW:5:20 HEADCROP:13 MINLEN:36

解釋一下這些參數(shù)

• trimmomatic PE: 運(yùn)行 Trimmomatic 的命令行接口，并指定輸入文件是配對的 PE (paired-end) reads

• -threads 1: 指定使用 1 個線程進(jìn)行處理

• -phred33: 指定輸入文件的 Phred 質(zhì)量值編碼是 Phred-33

• SRR25909836_1.fastq.gz 和 SRR25909836_2.fastq.gz: 輸入的配對 FASTQ 文件，SRR25909836_1.fastq.gz 是第一條 reads 的文件，SRR25909836_2.fastq.gz 是第二條 reads 的文件

• -summary oryza_sativa.summary: 生成一個名為 oryza_sativa.summary 的文件，其中包含了程序運(yùn)行的詳細(xì)統(tǒng)計信息

• -baseout trimmed: 輸出文件的基本名稱，后面會加上一些后綴來區(qū)分不同的輸出文件

• LEADING:3: 移除每條 reads 開頭的低質(zhì)量堿基，如果堿基質(zhì)量低于 3，就被去除

• TRAILING:3: 移除每條 reads 末尾的低質(zhì)量堿基，如果堿基質(zhì)量低于 3，就被去除

• SLIDINGWINDOW:5:20: 對每條 reads 進(jìn)行滑動窗口的質(zhì)量修剪，如果窗口內(nèi)的平均質(zhì)量小于 20，窗口內(nèi)的所有堿基就會被去除，窗口的大小為 5

• HEADCROP:13: 去除每條 reads 開頭的 13 個堿基

• MINLEN:36: 去除所有長度小于 36 的 reads

-phred33 之后的兩個是正向和反向的測序文件

會看到五個輸出文件 SRR25909836_1P.fastq.gz、SRR25909836_2P.fastq.gz、SRR25909836_1U.fastq.gz、SRR25909836_2U.fastq.gz、oryza_sativa.summary

• SRR25909836_1P.fastq.gz、SRR25909836_2P.fastq.gz 這兩個文件包含那些在兩個末端都通過質(zhì)量控制的序列配對，這兩個文件是過濾后的測序數(shù)據(jù)，可以用這兩個文件進(jìn)行后續(xù)的測序數(shù)據(jù)分析

• SRR25909836_1U.fastq.gz、SRR25909836_2U.fastq.gz 兩個文件包含那些只有一端通過質(zhì)量控制的序列，這就是我們說的"單一通過"（unpaired pass）輸出

• summary 提供了關(guān)于 Trimmomatic 運(yùn)行的詳細(xì)信息，它包括以下幾類信息

• 輸入的讀數(shù)量

• 裁剪因引物/接頭污染或低質(zhì)量得分而去除的讀數(shù)量

• 由于長度不足而被丟棄的讀物數(shù)量

• 輸出的讀數(shù)量

hisat2 比對

hisat2 可以快速準(zhǔn)確地將測序得到的 RNA 片段（reads）比對到參考基因組，從而確定這些 RNA 片段在基因組上的精確位置，進(jìn)一步可以用于基因表達(dá)量定量，剪接位點的檢測等多種 RNA-Seq 分析任務(wù)

建立索引

hisat2 需要一個 index 索引才能進(jìn)行比對，hisat2 提供了一些 index，但很少，只有人類、小鼠等基因組的，可以在下面的 ftp 地址中進(jìn)行下載

ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data

由于這里我做的是水稻的，所以就需要自己建立索引，使用的是水稻的參考基因組序列，前面已經(jīng)下載好了，使用下列命令建立索引

hisat2-build -p 4 oryza_sativa.fa oryza_sativa

• -p 4 是指占用 4 線程

• oryza_sativa.fa 即前面我重命名的參考基因組序列文件

• oryza_sativa 就是參考基因組序列文件，fa 前面的內(nèi)容

建立完成之后就可以看到系統(tǒng)中多了一個 oryza_sativa 文件夾，cd 進(jìn)去就可以看到 8 個以 ht2 為拓展名的文件，要將參考基因組序列文件也放到這個文件夾

常用參數(shù)

• -x：參考基因組索引文件的前綴

• -1：正向測序文件

• -2：反向測序文件

• -S：輸出文件名（sam 文件）

• -P：線程數(shù)

• -t：打印加載索引文件和對齊讀取所需的時間

• –dta-cufflinks：出來的結(jié)果更適合 cufflinks 處理 （主要用于基因表達(dá)量的計算和差異表達(dá)基因的尋找）

• –no-unal：不記錄沒比對上的 reads

• –un-conc：在輸出文件寫出與參考基因組不一致的 reads 對

進(jìn)行比對

我使用的命令如下

hisat2 -x oryza_sativa/oryza_sativa -p 5 -1 trimmed_1P -2 trimmed_2P -S oryza_sativa.sam

注意-x 后跟索引文件，不加拓展名，保證 ht2 文件和 fa 文件的文件名一致即可，這里由于前面過濾后的序列是沒有拓展名的，所以會提示 Warning: Unsupported file format，不影響結(jié)果

運(yùn)行完畢后便得到 sam 文件，還會輸出一段信息

z -S oryza_sativa.sam
21675765 reads; of these:
  21675765 (100.00%) were paired; of these:
    2281601 (10.53%) aligned concordantly 0 times
    18927559 (87.32%) aligned concordantly exactly 1 time
    466605 (2.15%) aligned concordantly >1 times
    ----
    2281601 pairs aligned concordantly 0 times; of these:
      248174 (10.88%) aligned discordantly 1 time
    ----
    2033427 pairs aligned 0 times concordantly or discordantly; of these:
      4066854 mates make up the pairs; of these:
        2387493 (58.71%) aligned 0 times
        1619885 (39.83%) aligned exactly 1 time
        59476 (1.46%) aligned >1 times
94.49% overall alignment rate

這些輸出記錄包含以下信息：

• 總共有 21675765 個讀取序列

• 所有讀取序列中 100.00%都成對存在（即是成對端 (read pairs) 序列）

• 成對端序列中 10.53%的序列沒有成功比對到基因組上

• 87.32%的序列只比對到了基因組上的一個位置

• 2.15%的序列比對到了基因組上的多個位置

• 對于沒有成功比對的成對端序列，有 10.88%者的序列不一致地（非正確配對的）比對到了基因組上一個位置

• 有的序列無法一致地或不一致地比對，這些序列占所有沒有成功比對的成對端序列的 2033427 對，它們一共包含 4066854 個“pairs”序列

• 在這些“pairs”序列中，

• 58.71%的序列沒有比對到任何地方

• 39.83%的序列比對到了基因組上的一個位置

• 1.46%的序列比對到了基因組上的多個位置

• 整體上的比對成功率為 94.49%

在 RNA-Seq 分析中，比對成功率是一個重要的質(zhì)量控制指標(biāo)， 94.49%的比對成功率表明，絕大部分讀取序列都能夠成功地比對到基因組上，這表示 RNA-Seq 實驗和測序質(zhì)量都相對較好

samtools 排序壓縮

常用命令

• samtools view：將 SAM 格式文件轉(zhuǎn)換為 BAM 格式，或者執(zhí)行過濾和查看操作

• -b（輸出 BAM 格式）

• -S（輸入為 SAM 格式）

• -h（輸出頭文件）

• samtools sort：對 BAM 文件進(jìn)行排序

• -o（輸出至文件）

• -n（按 read 名稱排序）

• samtools index：為排序后的 BAM 文件建立索引

• samtools faidx：為 fasta 格式的參考序列建立索引，并可以快速抽取序列

• samtools tview：文本模式下查看比對結(jié)果

• samtools flagstat：提供 BAM 文件的比對統(tǒng)計信息，如總 reads 數(shù)，映射的 reads 數(shù)等

• samtools mpileup：生成 mpileup 文件用于隨后的突變檢測

• -u（生成 BCF）

• -g（生成 BCF 或 VCF）

• -f（參考序列）

排序壓縮

我使用的是以下命令

samtools sort -n -@ 5 oryza_sativa.sam -o oryza_sativa

運(yùn)行完成后會得到一個 bam 文件

featureCounts 生成基因計數(shù)表

常用參數(shù)

• -a：注釋文件的路徑，格式為 GTF 或 GFF

• -o：輸出文件的路

• -p：如果序列數(shù)據(jù)來自于 paired-end RNA-seq，使用此參數(shù)

• -B：只保留同時映射到同一基因組特征（例如，基因或外顯子）的片段

• -C：不包括由基因間跨越引起的 reads

• -s：設(shè)置 strand-specific read 計數(shù)，如果使用 0，表示非特異性；如果使用 1，表示第一鏈特異性；如果使用 2，表示第二鏈特異性

• -Q：包含比對質(zhì)量達(dá)到特定閾值的 reads，默認(rèn)值為 10

• -g：制定 GTF 特性類型，例如，gene, transcript 等

• -F：制定輸入文件格式，例如，BAM, SAM 等

計數(shù)統(tǒng)計

我使用的是以下命令

featureCounts -T 5 -t exon -g Name -a oryza_sativa.gff3 -o gene_name.counts -p oryza_sativa

oryza_sativa.gff3 就是最初下載的注釋文件，如果要統(tǒng)計多個文件的話，在-p 后面跟上就可以，會生成 gene_name.counts、gene_name.counts.summary 兩個文件，

gene_name.counts.summary 文件是計數(shù)統(tǒng)計情況

gene_name.counts 文件是基因的具體信息

我這里只有一組數(shù)據(jù)，所以數(shù)量統(tǒng)計也只有一列，通常做 RNA-Seq 時是需要多組數(shù)據(jù)進(jìn)行分析的

參考

• RNASeq 原始數(shù)據(jù)質(zhì)量控制-FastQC

• 20160410 測序分析——使用 FastQC 做質(zhì)控

• RNA-seq 轉(zhuǎn)錄組數(shù)據(jù)分析丨一套完整的案例流程

• 生信軟件 | Trimmomatic （測序數(shù)據(jù)質(zhì)控）

• SRA 文件的下載（prefetch）和解壓 SRA 文件（fastq-dump）

• 從零開始的 RNASeq 教程（三）獲得基因計數(shù)表格

• 轉(zhuǎn)錄組差異表達(dá)分析–免比對工具 kallisto