對比軟件有很多，首先你可以收藏一下，因為我們是給你介紹的，請統(tǒng)一使用，但要了解它的用法。

直接去主頁下載index文件，然后對比fastq格式的reads，得到sam文件。

然后用它轉(zhuǎn)成bam文件，并排序（注意N排序和P排序的區(qū)別）并索引，加載到IGV中，截圖幾個基因看看！

順便對bam文件做一個簡單的QC，參考live 系列。

接下來的四篇文章基本都是規(guī)劃工作，從這篇文章開始我們才進入了RNA-Seq數(shù)據(jù)分析的核心部分。

比較還是不比較

在比較之前，我們先要明白比較的目的是什么？ RNA-Seq數(shù)據(jù)比較和DNA-Seq數(shù)據(jù)比較有什么區(qū)別？

RNA-Seq 數(shù)據(jù)分析有多種類型，例如尋找差異表達基因或?qū)ふ倚碌奶娲艚印?如果您正在尋找不同表達的基因，您只需要確定不同的讀取技術(shù)即可。 我們可以使用bwa這樣的比較工具，或者像這樣的免對齊工具，但前者更快。

如果您需要尋找新的或替代的 RNA 剪接，您需要 STAR 等工具來查找剪接位點。 由于RNA-Seq與DNA-Seq不同，當(dāng)DNA轉(zhuǎn)錄為mRNA時，內(nèi)含子被部分去除。 因此，如果mRNA倒轉(zhuǎn)的cDNA無法與參考序列進行比較，則會將其分離并重新比對，以確定中間是否有內(nèi)含子。

工具選擇

2016年-的綜述中提到，目前RNA數(shù)據(jù)分析有三種策略。 當(dāng)時的工具也主要使用，STAR和. 其中，其作者已推薦使用HISAT來代替。

最近發(fā)表了一篇題為-RNA-seq的文章——被譽為史上最完整的RNA-Seq數(shù)據(jù)分析流程，也是我一直想做的事情，但他們所做的卻超出了我的想象。 文章中根據(jù)參考基因組分析轉(zhuǎn)錄本所用的工具是、和STAR，推論是找到正確率最高的，總量小于和STAR。 從這里可以看出，第二類錯誤（接受錯誤）相對較小，而第一類錯誤（放棄正確）則較高。

就唯一的比較而言，STAR是兩者中最好的，主要是因為它不會像STAR和STAR那樣在PE匹配失敗時強行將SE與基因組進行比較。 而且在處理較長reads和較短reads的情況下，STAR的穩(wěn)定性也是最好的。

速度方面，比STAR和平均速度快2.5~100倍。

如果你是學(xué)習(xí)RNA-Seq數(shù)據(jù)分析，里??面提到的兩篇文檔一定要讀3遍以上，建議每隔一段時間復(fù)習(xí)一下。 而且就比較工具而言，基本上就選有STAR的。

下載索引

首先問自己一個問題，為什么比較的時候需要使用索引？ 這里強烈推薦大家閱讀Jimmy寫的算法原理探索。 而且這只是一個建議，你不需要真正讀它，因為你無論如何也無法理解它。

定量測序遇到的第一個問題是，如果在合理的時間內(nèi)將數(shù)千甚至數(shù)億的reads與參考基因組進行比對，并且保證錯誤率在可接受的范圍內(nèi)。 為了提高比對速度，需要通過BWT算法將參考基因組序列轉(zhuǎn)換為索引，而我們比對的序列可能是索引的子集。 事實上，轉(zhuǎn)錄組比較還考慮到了可變剪接的情況，所以比較復(fù)雜。

因此，我們不會直接將讀數(shù)發(fā)布到基因組中，而是將讀數(shù)與索引進行比較。 通常有現(xiàn)成的索引可供人類使用。 我建議你下載現(xiàn)成的。 之前我也嘗試過使用服務(wù)器自己創(chuàng)建索引，所花費的時間讓我懷疑自己的人生。

cd referece && mkdir index && cd index
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
tar -zxvf hg19.tar.gz

我認(rèn)為筆記本的配置相當(dāng)不錯，或者如果沒有現(xiàn)成的索引，可以通過以下方式創(chuàng)建

# 其實hisat2-buld在運行的時候也會自己尋找exons和splice_sites，但是先做的目的是為了提高運行效率
extract_exons.py gencode.v26lift37.annotation.sorted.gtf > hg19.exons.gtf &
extract_splice_sites.py gencode.v26lift37.annotation.gtf > hg19.splice_sites.gtf &
# 建立index， 必須選項是基因組所在文件路徑和輸出的前綴
hisat2-build --ss hg19.splice_sites.gtf --exon hg19.exons.gtf genome/hg19/hg19.fa hg19

我的是i7-7700處理器，顯存64G，運行的資源效率如下：

很快就會比較

基本用法是[]*-x{-1-2|-U}[-S]，基本上提供了索引、PE數(shù)據(jù)或SE數(shù)據(jù)存儲位置的位置。 但其他可選參數(shù)是進步的好方法。 菜鳥使用默認(rèn)參數(shù)。

由于RNA-Seq數(shù)據(jù)來自~，6個樣本的PE數(shù)據(jù)，即有6個循環(huán)，可以編寫腳本，或者直接在命令行運行

下面命令運行的目錄是/mnt/f/Data/，我的參考基因組的索引數(shù)據(jù)存放在/mnt/f/Data//index/hg19/，RNA-seq數(shù)據(jù)存放在/ mnt/f/Data/RNA-Seq 下降。 比對結(jié)果將存儲在/mnt/f/Data/RNA-Seq/中

mkdir -p RNA-Seq/aligned
for i in `seq 57 62`
do
 ? ?hisat2 -t -x reference/index/hg19/genome -1 RNA-Seq/SRR35899${i}_1.fastq.gz -2 SRR35899${i}_2.fastq.gz -S RNA-Seq/aligned/SRR35899${i}.sam &
done

& 會將任務(wù)丟到后臺，所以這3個比較程序會同時執(zhí)行，如果CPU和顯存無法承受，就將&一一去掉。 比較步驟消耗了大量的顯存資源，這是由于比較工具將索引數(shù)據(jù)加載到顯存中造成的。 我有64G顯存，理論上可以同時處理20個PE數(shù)據(jù)。 在我的筆記本電腦配置下，同時完成此步驟大約需要2個小時。

基本參數(shù)說明

對比數(shù)據(jù)時，可以悄悄讀取額外的選項，主要分為以下幾個部分

注：soft表示比對的reads中只有部分與參考序列匹配，有部分不匹配。 也就是說，一個100bp的read匹配下一個20bp或者前面的20bp，或者前面20bp的比較效果不是很好。

因此，影響比較結(jié)果的是比較選項、評分選項、可變剪切選項和PE選項。 我應(yīng)該在有生之年寫一篇文章來介紹這個選項對結(jié)果的影響。

輸出結(jié)果

比較后dnastar序列比對，將輸出以下結(jié)果。 分析顯示，所有數(shù)據(jù)都是100%，96.68%的配對數(shù)據(jù)沒有進行過一次比較，1.23%的數(shù)據(jù)是唯一比較，2.09%是多次比較。 之后，96.68%的數(shù)據(jù)沒有一次對比數(shù)據(jù)。 如果不按順序，則有0.05%的比較。 然后將其余數(shù)據(jù)與推拉數(shù)據(jù)進行比較，95.20%的數(shù)據(jù)沒有進行比較，3.60%的數(shù)據(jù)進行了一次比較，1.20%的數(shù)據(jù)進行了多次比較。 零和零的總和是對比的8%！ ！ ！

這個整體對比率讓我開始懷疑自己的人生。 這怎么可能？ 我查看了輸出記錄，發(fā)現(xiàn)幾個結(jié)果的比對率都超過了90%。 我想了想，這個實驗好像是人和老鼠都做的。 于是我又去GEO查了一下記錄，突然發(fā)現(xiàn)自己差點翻船了。

圖9-15是人293細胞(9-11)或小鼠ES細胞(12-15)中的mRNA序列。

同時我也反思了錯誤的原因。 我默認(rèn)這個實驗重復(fù) 3 次 KO 和非 KO。 雖然文章的實驗設(shè)計不是這樣的，但是可以看出理解實驗設(shè)計是非常重要的。 。

mkdir -p RNA-Seq/aligned
for i in `seq 56 58`
do
 ? ?hisat2 -t -x reference/index/hg19/genome -1 RNA-Seq/SRR35899${i}_1.fastq.gz -2 SRR35899${i}_2.fastq.gz -S RNA-Seq/SRR35899${i}.sam &
done

以上是修改后的代碼

三軸

SAM（/）數(shù)據(jù)格式是目前測序中存儲比對數(shù)據(jù)的標(biāo)準(zhǔn)格式。 事實上，它可以用來存儲不匹配的數(shù)據(jù)。那么，SAM的格式說明

目前處理SAM格式的工具主要是fúLi手寫的。 不僅有C語言版本，還有Java、Pysam、lisp cl-sam等版本。 主要功能如下：

最常用的三個軸是格式轉(zhuǎn)換、排序和索引。 進階教程是看文檔增強。

for i in `seq 56 58`
do
 ? ?samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam
 ? ?samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam
 ? ?samtools index SRR35899${i}_sorted.bam
done

筆記

Jimmy說，我們仔細判斷了sam排序的兩種形式的區(qū)別，于是我截取了上面的100行數(shù)據(jù)，分別排序，檢查結(jié)果。

head -1000 SRR3589957.sam > test.sam
samtools view -b ?test.sam > test.bam
samtools view test.bam | head

默認(rèn)排序

samtools sort test.bam default
samtools view default.bam | head

，orby 在使用時讀取名稱

samtools sort -n test.bam sort_left
samtools view sort_left.bam | head

另外說說默認(rèn)按照染色體位置排序，-n參數(shù)是按照read name排序。 其實還有一個-t可以按TAG排序。

談?wù)撚^點

三板夫的view是一個非常實用的子命令，除了前面的格式轉(zhuǎn)換之外，還可以進行數(shù)據(jù)的提取和提取。

例如，提取1號染色體1234區(qū)域的比對讀數(shù)

samtools view SRR3589957_sorted.bam chr1:1234-123456 | head

例如在搭配標(biāo)志（0.1.19版本沒有）中，使用-f或-F參數(shù)來提取不同匹配情況的讀取。

flag是描述reads比較的標(biāo)志。 有 12 種類型的標(biāo)志可以一起使用。

0x1 ? ?PAIRED ? ?paired-end (or multiple-segment) sequencing technology
0x2 ? ?PROPER_PAIR ? ?each segment properly aligned according to the aligner
0x4 ? ?UNMAP ? ?segment unmapped
0x8 ? ?MUNMAP ? ?next segment in the template unmapped
0x10 ? ?REVERSE ? ?SEQ is reverse complemented
0x20 ? ?MREVERSE ? ?SEQ of the next segment in the template is reverse complemented
0x40 ? ?READ1 ? ?the first segment in the template
0x80 ? ?READ2 ? ?the last segment in the template
0x100 ? ?SECONDARY ? ?secondary alignment
0x200 ? ?QCFAIL ? ?not passing quality controls
0x400 ? ?DUP ? ?PCR or optical duplicate
0x800 ? ?SUPPLEMENTARY ? ?supplementary alignment

你可以先看看整體情況

samtools flagstat SRR3589957_sorted.bam

也就是說，如果我想過濾剛剛配對的read，我需要使用0x10

samtools view -b -f 0x10 SRR3589957_sorted.bam chr1:1234-123456 ?> flag.bam
samtools flagstat flag.bam

我有生之年應(yīng)該寫一篇文章介紹一下。

比較質(zhì)量控制 (QC)

還是在-上，提到人類基因組的比對率應(yīng)該是70%到90%，并且多讀（multi-reads）的次數(shù)應(yīng)該少一些。 此外，對齊的外顯子和對齊的鏈（讀取一個外顯子和）的覆蓋范圍應(yīng)該一致。

常用的工具有

我們來使用RSeQC，雖然使用的是寫入工具，但是自然友好dnastar序列比對，安裝也很方便。

# Python2.7環(huán)境下
pip install RSeQC

總共有以下幾個文件，根據(jù)名稱就可以知道功能是什么。

首先對bam文件進行統(tǒng)計分析。 從結(jié)果來看，符合70~90的對比率要求。

bam_stat.py -i SRR3589956_sorted.bam

基因組覆蓋度的QC需要提供bed文件，可以直接從RSeQC網(wǎng)站下載，也可以用gtf轉(zhuǎn)換

read_distribution.py -i RNA-Seq/aligned/SRR3589956_sorted.bam -r reference/hg19_RefSeq.bed

IGV視圖

加載參考序列、注釋和 BAM 文件，然后查看。

大家好奇為什么一開始是五嗎，下面會貼出答案，原文是我的地址。

()

如有侵權(quán)請聯(lián)系刪除！