有很多比較軟件。 首先，請收集它們。 因為我們是來向您介紹它們的，所以請統(tǒng)一使用它們并了解如何使用它們。

只要去首頁下載index文件，然后對比fastq格式的reads就可以得到sam文件了。

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

順便對bam文件做一下簡單的QC，參考我的基因組系列直播。

前四篇文章基本上都是準備工作。 從本文開始，我們進入RNA-Seq數(shù)據(jù)分析的核心部分。

比較還是不比較

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

RNA-Seq數(shù)據(jù)分析可以分為多種類型，例如尋找差異表達基因或?qū)ふ倚碌倪x擇性剪接。 如果我們要尋找差異表達的基因，我們只需要確定不同的read技術(shù)即可。 我們可以使用bwa等比較工具，或者align-free工具，后者速度更快。

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

工具選擇

在 2016 年的一篇評論 A of best for RNA-seq data 中，提到目前 RNA 數(shù)據(jù)分析有三種策略。 當(dāng)時主要使用的工具是STAR和. 其中，其作者已推薦使用HISAT來代替。

近日，一篇題為《by a Broad-RNA-seq》的文章發(fā)表——被稱為史上最全面的RNA-Seq數(shù)據(jù)分析流程。 這也是我一直想做的，但他們做得很好。 超出我的想象。 文章中基于參考基因組進行轉(zhuǎn)錄本分析所使用的工具是、和STAR。 結(jié)論是準確率最高，但總數(shù)少于STAR。 從這里可以看出，第2類錯誤（接受謊言）相對較少，但第1類錯誤（拋棄真相）數(shù)量較多。

就獨特的比較而言，STAR是三者中最好的，主要是因為它在PE比較失敗時不會將SE強制到基因組上，就像和一樣。 而且，在處理長讀和短讀的不同情況時，STAR的穩(wěn)定性也是最好的。

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

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

下載索引

首先，問自己一個問題，為什么比較時需要使用索引？ 強烈推薦大家閱讀Jimmy寫的算法原理探索和算法原理探索。 但這只是一個建議，你不需要實際看到，反正你也不會理解。

高通量測序遇到的第一個問題是，必須在合理的時間內(nèi)將數(shù)千甚至數(shù)億的reads映射到參考基因組上，并且錯誤率必須在可接受的范圍內(nèi)。 為了提高比對速度，需要通過BWT算法將參考基因組序列轉(zhuǎn)換為索引，而我們比對的序列實際上是索引的子集。 當(dāng)然，轉(zhuǎn)錄組比較還需要考慮選擇性剪接，因此比較復(fù)雜。

因此，我們不是將讀數(shù)直接粘貼到基因組上，而是將讀數(shù)與索引進行比較。 人類索引一般都是現(xiàn)成的。 我建議您下載現(xiàn)成的。 我曾經(jīng)嘗試使用服務(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

如果你覺得你的電腦配置還可以，或者沒有現(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處理器，內(nèi)存64G。 運行資源效率如下：

形式比較

基本用法是 []* -x {-1 -2 | -U } [-S]，基本上提供了索引、PE數(shù)據(jù)或SE數(shù)據(jù)存儲位置的位置。 然而，其他可選參數(shù)是一個很好的進步方式。 如果您是新手，則使用默認參數(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 測序下。 比較結(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ù)扔到后臺，因此這三個比較程序會同時執(zhí)行。 如果CPU和內(nèi)存無法承受，請將&一一刪除。 比較步驟消耗大量內(nèi)存資源。 這是由于比較工具必須將索引數(shù)據(jù)放入內(nèi)存而導(dǎo)致的。 我有64G內(nèi)存，理論上可以同時處理20個PE數(shù)據(jù)。 在我的電腦配置下dnastar序列比對，同時完成這一步大約需要2個小時。

基本參數(shù)說明

在數(shù)據(jù)比較過程中，您可以有額外的安靜讀取選項，主要分為以下幾個部分：

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

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

輸出結(jié)果

比較后，將輸出以下結(jié)果。 解讀一下，所有數(shù)據(jù)都是100%，96.68%的配對數(shù)據(jù)不是一次性比較的，1.23%的數(shù)據(jù)比較是唯一比較，2.09%是多重比較。 那么96.68%的人根本沒有數(shù)據(jù)可以比較。 如果不按順序，就會有0.05%的比較。 如果再用單端數(shù)據(jù)進行比較，則95.20%的數(shù)據(jù)沒有進行比較，3.60%的數(shù)據(jù)進行了一次比較，1.20%的數(shù)據(jù)進行了多次比較。 零和零的總和加起來就是8%的對比！ ！ ！

這個整體對比率讓我對生活產(chǎn)生了質(zhì)疑。 怎么可能？ 我查看了輸出記錄，發(fā)現(xiàn)有幾個結(jié)果的比對率都在90%以上。 我想了想，發(fā)現(xiàn)這個實驗似乎是在人類和老鼠身上都做過的。 于是又去GEO查了一下記錄，突然發(fā)現(xiàn)自己差點翻車了。

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

同時我也反思了錯誤的原因。 我默認 KO 和非 KO 各重復(fù) 3 次。 事實上，文章的實驗設(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ù)的標準格式。 當(dāng)然也可以用來存儲未對齊的數(shù)據(jù)。因此，SAM的格式說明

目前處理SAM格式的工具主要是Heng 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

默認排序

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

按 排序，或使用 -n 時按讀取名稱排序

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

也就是說，默認是根據(jù)染色體位置排序dnastar序列比對，而-n參數(shù)則是根據(jù)read name排序。 當(dāng)然還有一個-t根據(jù)TAG排序。

談?wù)効捶?/p>

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

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

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

例如，使用flag（0.1.19版本中不可用）和 ，使用-f或-F參數(shù)來提取不同匹配情況的reads。

Flag是一種描述讀取比較情況的標記。 共有12種，可以組合使用。

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

換句話說，如果我想過濾完全匹配的reads，我需要使用0x10

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

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

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

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

常用的工具包括

讓我們使用 RSeQC。 畢竟使用自己寫的工具本來就很友好，安裝也很方便。

# 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)系刪除！