superqun原創(chuàng)

一决记、流程概括

• RNA-seq的原始數(shù)據(jù)（raw data）的質(zhì)量評估
• linux環(huán)境和R語言環(huán)境
• raw data的過濾和清除不可信數(shù)據(jù)（clean reads）
• reads回帖基因組和轉(zhuǎn)錄組（alignment）
• 計(jì)數(shù)（count ）
• 基因差異分析（Gene DE）
• 數(shù)據(jù)的下游分析

二、準(zhǔn)備工作

• 學(xué)習(xí)illumina公司測序原理
• 測序得到的fastq文件
• 注釋文件和基因組文件的準(zhǔn)備

1. fastq測序文件

在illumina的測序文件中河劝，采用雙端測序（paired-end），一個樣本得到的是seq_1.fastq.gz和seq_2.fastq.gz兩個文件，每個文件存放一段測序文件橄唬。在illumina的測序的cDNA短鏈被修飾為以下形式（圖源見水印）：

inllumina公司測序的cDNA短鏈形式

兩端的序列是保護(hù)堿基（terminal sequence）参歹、接頭序列（adapter）仰楚、索引序列（index）、引物結(jié)合位點(diǎn)（Primer Binding Site）：其中 adapter是和flowcell上的接頭互補(bǔ)配對結(jié)合的；index是一段特異序列僧界，加入index是為了提高illumina測序儀的使用率侨嘀，因?yàn)橥粋€泳道可能會測序多個樣品，樣品間的區(qū)分就是通過index區(qū)分捂襟。參考：illumina 雙端測序（pair end）咬腕、雙端測序中read1和read2的關(guān)系。

在illumina公司測得的序列文件經(jīng)過處理以fastq文件協(xié)議存儲為*.fastq格式文件葬荷。在fastq文件中每4行存儲一個read涨共。
第一行：以@開頭接ReadID和其他信息，分別介紹了
第二行：read測序信息
第三行：規(guī)定必須以“+”開頭宠漩，后面跟著可選的ID標(biāo)識符和可選的描述內(nèi)容煞赢，如果“+”后面有內(nèi)容，該內(nèi)容必須與第一行“@”后的內(nèi)容相同
第四行：每個堿基的質(zhì)量得分哄孤。記分方法是利用ERROR P經(jīng)過對數(shù)和運(yùn)算分為40個級別分別與ASCII碼的第33號!和第73號I對應(yīng)照筑。用ASCII碼表示堿基質(zhì)量是為了減少文件空間占據(jù)和防止移碼導(dǎo)致的數(shù)據(jù)損失。fastq文件預(yù)覽如下：

在這里插入圖片描述

2.注釋文件和基因組文件的獲取

• 基因組獲取方式：可以從NCBI瘦陈、NCSC凝危、Ensembl網(wǎng)站或者檢索關(guān)鍵詞“hg38 ftp UCSC” 人類基因組hg38.fa.gz大概是938MB左右。文件獲取可以點(diǎn)擊網(wǎng)站下載晨逝。可以通過云盤的離線下載來加速下載進(jìn)程
• 基因組的選擇：以Ensembl網(wǎng)站提供的基因組為例蛾默，比對用基因組應(yīng)該選擇Homo_sapiens.GRCh38.dna.primary_assembly.fa
• Ensembl基因組的不同版本詳見README和高通量測序數(shù)據(jù)處理學(xué)習(xí)記錄（零）：NGS分析如何選擇合適的參考基因組和注釋文件

三、軟件安裝

• 安裝方式：軟件安裝可以通過例如apt-get捉貌、miniconda等方式來安裝支鸡。由于miniconda的便捷行，使用conda進(jìn)行如下軟件的安裝趁窃。
• 軟件列舉
  質(zhì)控：fastqc ,multiqc , trimmomatic, cutadapt, trim-galore
  比對：star , hisat2 , tophat , bowtie2 , bwa , subread
  計(jì)數(shù)：htseq , bedtools, salmon
• miniconda的安裝：

1. 可以通過點(diǎn)擊清華大學(xué)開源軟件站或者檢索“清華大學(xué) conda”訪問鏡像網(wǎng)站（清華鏡像站因?yàn)榉?wù)器在中國訪問速度比較快）牧挣，點(diǎn)擊Anoconda界面，選擇Miniconda下載安裝醒陆，windows在安裝好需要設(shè)置環(huán)境變量瀑构。（該鏡像將于近日停用可以更換為中科大或者其他源）

清華大學(xué)宣布停止Anaconda鏡像服務(wù)

2. linux測試Miniconda的安裝：conda -v
3. 創(chuàng)建名為rna的環(huán)境變量：conda create -n rna python=2（許多軟件依賴python2環(huán)境）
4. 配置conda，添加鏡像源頭：輸入如下命令(暫時可用刨摩，但是該鏡像將于近日停用寺晌。屆時會更新新的鏡像地址和設(shè)置方法)

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
conda config --set show_channel_urls yes
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/

• 軟件安裝：conda install <software>會自動安裝軟件和軟件環(huán)境。值得注意的是需要在rna的環(huán)境變量下安裝以上軟件澡刹。激活rna環(huán)境變量的代碼：

source activate rna

四呻征、質(zhì)量匯報生成與讀取

1.fastq質(zhì)量匯報

使用命令fastqc -o <output dir> <seqfile1,seqfile2..>來進(jìn)行質(zhì)量報告。每個fastqc文件會獲得一個質(zhì)量分析報告罢浇，來描述此次RNA-seq的測序質(zhì)量陆赋。
獲取質(zhì)量報告如圖：

在這里插入圖片描述

Basic Statistics

從read水平來總覽边篮，判斷測序質(zhì)量。
Encoding ：測序平臺的版本奏甫，因?yàn)椴煌姹镜?error p的計(jì)算方法不一樣戈轿。
Total sequence：測序深度。一共測序的read數(shù)阵子。是質(zhì)量分析的主要參數(shù)思杯。
Sequence length：測序長度。
%GC：GC堿基含量比挠进，一般是物種特異性色乾，比如人類是42%左右。

Perbase sequence quality

橫坐標(biāo)： 第1-100個測序得到的堿基
縱坐標(biāo)： 測序質(zhì)量評估领突。這里的Q=-10*lg10(error P),即20%代表1%的錯誤讀取率暖璧，30%代表0.1%的錯誤讀取率
箱型圖： 紅色線，是某個順序下測序堿基所有測序質(zhì)量的中位數(shù)君旦。黃色塊澎办，是測序質(zhì)量在25%-75%區(qū)域。藍(lán)色線金砍，平均數(shù)局蚀。
一般要求： 測序箱型圖10%的線大于Q=20。Q20過濾法恕稠。

per tail sequence quality

橫坐標(biāo)：同上琅绅。
縱坐標(biāo)：tail的index編號。
目的：防止測序過程中某些tail受不可控因素測序質(zhì)量低鹅巍。
標(biāo)準(zhǔn)：藍(lán)色表示質(zhì)量高千扶，淺色或暖色表示質(zhì)量低，后續(xù)的分析可以去除低質(zhì)量tail骆捧。

在這里插入圖片描述

Per sequence quality scores

從read的總體測序質(zhì)量分布來判定此次的測序質(zhì)量澎羞，是質(zhì)量分析的重要標(biāo)準(zhǔn)之一。
橫坐標(biāo)：表示read的測序質(zhì)量Q=-10*lg10(error P)凑懂。
縱坐標(biāo)：表示在該Q值下的read 的數(shù)量
標(biāo)準(zhǔn)：需要集中在高分區(qū)

在這里插入圖片描述

Per base sequence content

橫坐標(biāo)：1-100的測序堿基位置
縱坐標(biāo)：堿基百分比
標(biāo)準(zhǔn)：理論上煤痕，ATCG堿基的分布應(yīng)該差別不大，即四條線應(yīng)該大致平行狀態(tài)接谨。如果AT或CG差異超過10%，此項(xiàng)檢測是危險的塘匣。一般是測序機(jī)器前幾個堿基測序時候因?yàn)闋顟B(tài)調(diào)整導(dǎo)致測序略有偏差脓豪，如果前幾個堿基偏差較大，可以在后期將前幾個堿基切掉忌卤。

在這里插入圖片描述

Sequence Length Distribution

統(tǒng)計(jì)read的堿基長度扫夜，本例理論上測序應(yīng)該全是100bp。
橫坐標(biāo)：是read的堿基長度
縱坐標(biāo)：是該長度下的read數(shù)量

在這里插入圖片描述

Per sequence GC content

橫坐標(biāo)：每個read的平局GC含量占比
縱坐標(biāo)：一定GC比下的read數(shù)
標(biāo)準(zhǔn)：藍(lán)色是理論值，紅色是真實(shí)值笤闯。兩者接近是比較好的狀態(tài)堕阔。如果有雙峰，可能混有了其他物種的DNA序列颗味。

在這里插入圖片描述

Adapter Content

一般測序在初步生成fastq文件時候超陆，adapter會被去除，但是有的會沒有去除或者遺漏部分adapter浦马。所以這一步是檢測RNA-seq測序過程中adapter是否去除时呀。如果沒有去除會嚴(yán)重影響后續(xù)的比對工作。沒有去除的adapter在質(zhì)量處理環(huán)節(jié)會被處理掉晶默。

在這里插入圖片描述

2. multiqc質(zhì)量報告

multiqc可以對幾個fastqc報告文件進(jìn)行總結(jié)并匯總到一個報告文件中谨娜，以更直觀到防止展示。使用方法

multiqc <analysis directory> 

在這里插入圖片描述

五磺陡、數(shù)據(jù)處理

數(shù)據(jù)處理內(nèi)容：fiter the bad quality reads and remove adaptors.
處理軟件：數(shù)據(jù)到處理可以使用多款軟件趴梢，trim_galore在各文獻(xiàn)中表現(xiàn)良好。

1.trim_galore 的使用方法

trim_galore：可以處理illumina币他，nextera3垢油，smallRNA測序平臺的雙端和單端數(shù)據(jù)，包括去除adapter和低質(zhì)量reads圆丹。
trim_galore的參數(shù)： trim_galore的參數(shù)在處理過程比較重要：

trim_galore [options] <filename>
--quality<int>  #設(shè)定phred quality閾值滩愁。默認(rèn)20（99%的read質(zhì)量），如果測序深度較深辫封，可以設(shè)定25
--phred33       #設(shè)定記分方式硝枉，代表Q+33=ASCII碼的方式來記分方式。這是默認(rèn)值倦微。
--paired          # 對于雙端結(jié)果妻味，一對reads中若一個read因?yàn)橘|(zhì)量或其他原因被拋棄，則對應(yīng)的另一個read也拋棄欣福。
--output_dir   #輸出目錄责球，需確保路徑存在并可以訪問
--length        #設(shè)定長度閾值，小于此長度會被拋棄拓劝。這里測序長度是100我設(shè)定來75雏逾，感覺有點(diǎn)浪費(fèi)
--strency     #設(shè)定可以忍受的前后adapter重疊的堿基數(shù)，默認(rèn)是1郑临。不是很明白這個參數(shù)的意義
-e<ERROR rate>  #設(shè)定默認(rèn)質(zhì)量控制數(shù)栖博，默認(rèn)是0.1，即ERROR rate大于10%的read 會被舍棄厢洞，如果添加來--paired參數(shù)則會舍棄一對reads
<filename>  #如果是采用illumina雙端測序的測序文件仇让，應(yīng)該同時輸入兩個文件典奉。

構(gòu)建命令：

trim_galore -output_dir clean --paired --length 75 --quality 25 --stringency 5 seq_1.fasq.gz seq_2.fastq.gz

處理需要花上一定時間和磁盤空間。得到處理后數(shù)據(jù)

在這里插入圖片描述

2. 整理后數(shù)據(jù)的質(zhì)量分析丧叽。

對過濾后對文件進(jìn)行質(zhì)量分析卫玖。觀察過濾結(jié)果。同樣使用fastqc和multiqc兩個軟件進(jìn)行質(zhì)量分析踊淳。得到結(jié)果如下：

在這里插入圖片描述

ENCFF108UVC_val_1_fastqc的質(zhì)量報告

觀察到總read數(shù)減小和總體read的質(zhì)量變高假瞬，小部分adapter也被去除。更具體過濾和trim_galore的數(shù)據(jù)處理情況可以在seq_trimming_report.txt中查看嚣崭。

SUMMARISING RUN PARAMETERS
==========================
Input filename: ENCFF108UVC.fastq.gz
Trimming mode: paired-end
Trim Galore version: 0.5.0
Cutadapt version: 1.18
Quality Phred score cutoff: 25
Quality encoding type selected: ASCII+33
Adapter sequence: 'AGATCGGAAGAGC' (Illumina TruSeq, Sanger iPCR; auto-detected)
Maximum trimming error rate: 0.1 (default)
Minimum required adapter overlap (stringency): 5 bp
Minimum required sequence length for both reads before a sequence pair gets removed: 75 bp
Output file will be GZIP compressed


This is cutadapt 1.18 with Python 2.7.6
Command line parameters: -f fastq -e 0.1 -q 25 -O 5 -a AGATCGGAAGAGC ENCFF108UVC.fastq.gz
Processing reads on 1 core in single-end mode ...
Finished in 1038.93 s (40 us/read; 1.50 M reads/minute).
=== Summary ===

Total reads processed:              26,038,229
Reads with adapters:                   714,205 (2.7%)
Reads written (passing filters):    26,038,229 (100.0%)

Total basepairs processed: 2,603,822,900 bp
Quality-trimmed:              82,577,636 bp (3.2%)
Total written (filtered):  2,513,138,030 bp (96.5%)

由報告可以知道處理的具體詳情笨触。

六、比對回帖

概況：使用處理后的fastq文件和基因組與轉(zhuǎn)錄組比對雹舀，確定在轉(zhuǎn)錄組或者基因組中的關(guān)系芦劣。在轉(zhuǎn)錄組和基因組的比對采取的方案不同。分別是ungapped alignment to transcriptome和Gapped aligenment to genome说榆。
軟件：hisat2和STAR在比對回帖上都有比較好的表現(xiàn)虚吟。有文獻(xiàn)顯示，hisat2在納偽較少但是棄真較多签财，但是速度比較快串慰。STAR就比對而言綜合質(zhì)量比較好，在長短reads回帖上都有良好發(fā)揮唱蒸。由于hisat2的速度優(yōu)勢邦鲫，選擇hisat2作為本次比對的軟件。
在比對之前首先要先進(jìn)行索引文件的獲取或者制作神汹。

1. 索引文件的獲取

• 不同的比對軟件構(gòu)建索引方式不同庆捺，所用的索引也不盡相同
• 索引文件可以去網(wǎng)站下載也可以自己構(gòu)建。但是索引構(gòu)建會比較費(fèi)時間屁魏。建立索引文件需要大約一個小時（MAC: 2.6 GHz Intel Core i5/ 8 GB 1600 MHz DDR3） 滔以。
• 網(wǎng)站下載hisat2基因組索引：http://ccb.jhu.edu/software/hisat2/index.shtml
• 本地索引文件構(gòu)建參考了CSDN@ Richard_Jolin的構(gòu)建過程

• 索引文件的格式如下，是由多個文件構(gòu)成氓拼，要保證索引文件的格式和名稱部分一致你画。

  
  
  在這里插入圖片描述

2. hisat2的比對回帖

使用hisat2回帖

公式構(gòu)建根據(jù)hisat2 的使用說明書構(gòu)建了以下公式：

hisat2  -p 6    -x <dir of index of genome>  -1  seq_val_1.fq.gz   -2  seq_val_2.fq.gz  -S  tem.hisat2.sam

參數(shù)說明：

-p #多線程數(shù)
-x #參考基因組索引文件目錄和前綴
-1 #雙端測序中一端測序文件
-2 #同上
-S #輸出的sam文件

說明：在比對過程中，hisat會自動將雙端測序匹配同一reads并在基因組中比對桃漾，最后兩個雙端測序生成一個sam文件坏匪。比對回帖過程需要消耗大量時間和電腦運(yùn)行速度和硬盤存儲空間。5G左右fastq文件比對回帖過程消耗大概一個小時呈队，生成了17G的sam格式文件剥槐。回帖完成會生成一個回帖報告宪摧。

在這里插入圖片描述

samtools 軟件進(jìn)行格式轉(zhuǎn)換

SAM文件和BAM文件
samtools 是針對比對回帖的結(jié)果——sam和bam格式文件的進(jìn)一步分析使用的軟件。sam格式文件由于體量過大，一般都是使用bam文件來進(jìn)行存儲击狮。由于bam文件是二進(jìn)制存儲所以文件大小比sam格式文件小許多筹陵，大約是sam格式體積的1/6 。
samtools將sam轉(zhuǎn)換bam文件

samtools view -S seq.sam -b > seq.bam  #文件格式轉(zhuǎn)換
samtools sort seq.bam -0 seq_sorted.bam  ##將bam文件排序
samtools index seq_sorted.bam  #對排序后對bam文件索引生成bai格式文件沿彭，用于快速隨機(jī)處理朽砰。

至此一個回帖到基因組對RNA-seq文件構(gòu)建完成。這個seq_sourted.bam文件可以通過samtools或者IGV( Integrative Genomics Viewer)獨(dú)立軟件進(jìn)行查看喉刘。在IGV軟件中載入seq_sourted.bam文件瞧柔。
可以很直觀清晰地觀察到reads在基因組中的回帖情況和外顯子與內(nèi)含子的關(guān)系。

在這里插入圖片描述

3.對回帖bam文件進(jìn)行質(zhì)量評估睦裳。

**samtools falgstate **：統(tǒng)計(jì)bam文件中比對flag信息造锅，然后輸出比對結(jié)果。
公式：

samtools flagstate seq_sorted.bam > seq_sorted.flagstate

結(jié)果如下

47335812 + 0 in total (QC-passed reads + QC-failed reads)
3734708 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
46714923 + 0 mapped (98.69% : N/A)
43601104 + 0 paired in sequencing

21800552 + 0 read1
21800552 + 0 read2
42216752 + 0 properly paired (96.82% : N/A)
42879780 + 0 with itself and mate mapped
100435 + 0 singletons (0.23% : N/A)
337412 + 0 with mate mapped to a different chr
308168 + 0 with mate mapped to a different chr (mapQ>=5)

七廉邑、count

計(jì)算RNA-seq測序reads對在基因組中對比對深度哥蔚。
計(jì)數(shù)工具：feature counts
公式構(gòu)建：

feature counts -T 6 -t exon -g gene_id -a <gencode.gtf> -o seq_featurecount.txt <seq.bam>

參數(shù)：

-g # 注釋文件中提取對Meta-feature 默認(rèn)是gene_id
-t # 提取注釋文件中的Meta-feature 默認(rèn)是 exon
-p #參數(shù)是針對paired-end 數(shù)據(jù)
-a #輸入GTF/GFF 注釋文件
-o #輸出文件

接下來是表達(dá)矩陣構(gòu)建。在R語言環(huán)境下分析蛛蒙。

</br>

共勉糙箍！歡迎大家踴躍交流，討論牵祟，質(zhì)疑深夯，批評。留言必回诺苹。

我想建立并管理一個高質(zhì)量的生信&統(tǒng)計(jì)相關(guān)的微信討論群咕晋，如果你想?yún)⑴c討論，可以添加微信：veryqun 筝尾。我會拉你進(jìn)群捡需，當(dāng)然有問題也可以微信咨詢我。