1. 程式人生 > >一次rna-seq的過程-知乎live轉

一次rna-seq的過程-知乎live轉

生成 trim 顯示 百分比 sheng id號 通用 cat clas

數據分析流程

來自知乎孟浩巍的“快速入門生物信息學的”Live,超棒的~

技術分享圖片

首先是質控部分,使用fastqc進行對結果分析。

對於Illumia二代測序的結果質控包括兩個方面,去掉測序質量不好的序列,即Quality Control;二是需要去掉連在玻璃上的短的接頭,cut adaptor。

技術分享圖片

-t 8表示調用8個核心去運算。

之後,對每一個序列文件都生成一個zip和一個html文件。

例如:

技術分享圖片

那麽這2500000肯定是不同的基因,只不過這個機器的測序長度是150,所以所有的基因長度都是相同的。

技術分享圖片

對250萬reads,每個位點的Q做一個箱線圖,要求箱線值最低點高於20%,否則需要將那部分切除。

技術分享圖片

在145左右的的序列Q值較低測序不穩,所以這樣的序列145之後的全不要了。

技術分享圖片

這個是GC含量圖,通常A和T相同,C和G相同,但是前10bp不穩定,需要切除。

技術分享圖片

表示測序過程中測到的段序列的含量,橫軸是1-150bp,縱軸是百分比。由於某種原因導致測的絕大多數都是測的adaptor。

技術分享圖片

這個主要是衡量建庫水平,建庫中通常有6-8輪的PCR,但有時會出現過P的現象,當duplication過高的情況下,需要去dup。但是在RNA-seq裏通常是不去dup的。

②接下來,使用fastx_trimmer去頭去尾

技術分享圖片

zcat $fastq_1 | fastx_trimmer -f 11 -l 140 -z -o $out_fastq_1 &

zcat解壓縮,$fastq_1是輸入的第一個文件,這個文件解壓縮之後的結果給fastx_trimmer這個命令,

這個命令的參數-f是指first即保留的第一個bp(這裏前10bp剪切掉了);last即保留的最後一個bp(保留到第140bp),-z是壓縮命令,-o是輸出到這個文件裏。

//其中$:在bash裏表示當前是普通用戶;是變量引用操作符。a=10; echo $a會輸出10。

③使用cutadaptor去掉兩端的adaptor。

trimmer之後有一個去adaptor的過程,使用cutadaptor的軟件,

技術分享圖片

nohup cutadapt --times 1 -e 0.1 -0 3 
--quality-cutoff 6 -m 50 -a AGATCGGAAGAGC 
-A AGATCGGAAGAGC -o $out_fastq_1 
-p $out_fastq_2 $fastq_1 $fastq_2  >  $log_file 2>$1 &

//其中nohup:不掛斷地運行命令。

//2>$1:$1是傳遞給shell腳本的第一個參數;(轉自:https://www.cnblogs.com/kaituorensheng/p/4002697.html)

$# 是傳給腳本的參數個數
$0 是腳本本身的名字
$1 是傳遞給該shell腳本的第一個參數
$2 是傳遞給該shell腳本的第二個參數
$@ 是傳給腳本的所有參數的列表
$* 是以一個單字符串顯示所有向腳本傳遞的參數,與位置變量不同,參數可超過9個
$$ 是腳本運行的當前進程ID號
$? 是顯示最後命令的退出狀態,0表示沒有錯誤,其他表示有錯誤

//times 1一條序列只去一次Adaptor;-e 0.1在匹配時可以有10%的錯誤率;-O 3 adaptor序列必須和測序序列有3個堿基以上的overlap才可以;常用6;-m 50如果處理之後低於50的話就扔掉序列,短序列測序質量可能不是很好;-a和-A是Illumina常用的通用引物,之所以輸入兩個,是因為我是一個雙端測序的結果,需要對兩個文件內容進行分別去除,-a對應Reads1,-A對應reads2,$fasrq_1和_2是上一步的輸出;>最後是寫入log文件

一次rna-seq的過程-知乎live轉