需要用安裝好的sratoolkit把sra文件轉(zhuǎn)換為fastq格式的測序文件，并且用fastqc軟件測試測序文件的質(zhì)量！
作業(yè)，理解測序reads，GC含量，質(zhì)量值，接頭，index，fastqc的全部報告，搜索中文教程，并發(fā)在論壇上面。

數(shù)據(jù)解壓

之前下載了所有的數(shù)據(jù)，但只有樣本9_{15才是mRNA-Seq測序結(jié)果，其中9-11為人類293個AKAP9敲除細胞，12-15為小鼠的AKAP9敲除細胞。也就是只要解壓9}15就行了，即是SRR3589956 ~ SRR3589962.

先看一篇文章做過1000遍RNA-Seq的老司機告訴你如何翻車， 看一下事故現(xiàn)場，避開一些坑。
可以先用fastq-dump -h看一下幫助文件，分為如下幾個部分：

• 輸入： -A|--accession 序列號
• 處理中: Read Splitting, Full Spot Filters, Common Filters, Filters based on alignments, Filters for individual reads。 基本都是些過濾參數(shù)。不太常用
• 輸出： -O|--outdir 輸出文件夾， -Z|--stdout 輸出到標準輸出， --gzip/--bzip2 輸出為壓縮格式
• 多文件選項： 常用的就是--split-3
• 格式化： 分為序列， 質(zhì)量等，不常用

所以基本上命令即使如下用法, 如果你覺得自己的空間夠大就不需要用到--gzip參數(shù)。

for id in `seq 56 62`
do
    fastq-dump --gzip --split-3 -O /mnt/f/Data/RNA-Seq -A SRR35899${id}
done

先用 man seq了解下seq的用途，現(xiàn)在推薦用fasterq-dump進行解壓，速度更快

壓縮的文件不要著急解壓，有很多bash命令能夠直接用于壓縮文件，如zgrep,zcat,zless,zdiff等。

zcat SRR3589956_1.fastq.gz | head -n 4
@SRR3589956.1 D5VG2KN1:224:C4VAYACXX:5:1101:1159:2173 length=51
GGCGAGTGTAGGGCTGGCGCTGCCGGACGCGGTGCTAGTCGCCGGATGAAG
+SRR3589956.1 D5VG2KN1:224:C4VAYACXX:5:1101:1159:2173 length=51
B<BFBFBF0BFFFBFFBBFFIF<FFI<7<<BF<FFFFFFBB<BBBBBBBBB

用這些z-tools能夠節(jié)省大量磁盤空間。

QC basic concept

高通量測序之所以能夠能夠達到如此高的通量的原因就是他把原來幾十M，幾百M，甚至幾個G的基因組通過物理或化學的方式打算成幾百bp的短序列，然后同時測序。

因為測序過程是邊合成邊測序（SBS），所以在建庫的時候，短序列兩段會加一些固定的堿基用于橋式PCR擴增，這些固定的堿基就是adapter(接頭)。一般而言，還可以在接頭加一些tag（index），用于標識這個read來自于哪個物種。目前的單細胞測序為了省錢，譬如10X genomic技術(shù)，都是在一個pool里面加多種接頭。

在測序過程中，機器會對每次讀取的結(jié)果賦予一個值，用于表明它有多大把握結(jié)果是對的。從理論上都是前面質(zhì)量好，后面質(zhì)量差。并且在某些GC比例高的區(qū)域，測序質(zhì)量會大幅度降低。

因此，我們在正式的數(shù)據(jù)分析之前需要對分析結(jié)果進行質(zhì)控。Jimmy大神就發(fā)帖專門指出”要充分了解你的測序數(shù)據(jù)--論QC的重要性“，http://www.biotrainee.com/thread-324-1-1.html 。

Fastq文件格式說明

FASTQ文件每個序列通常為4行，分別為：

• Line 1 begins with a '@' character and is followed by a sequence identifier and an optional description (like a FASTA title line).
• Line 2 is the raw sequence letters.
• Line 3 begins with a '+' character and is optionally followed by the same sequence identifier (and any description) again.
• Line 4 encodes the quality values for the sequence in Line 2, and must contain the same number of symbols as letters in the sequence.

FASTQ的文件示例：

@DJB775P1:248:D0MDGACXX:7:1202:12362:49613 1:Y:18:ATCACG
TGCTTACTCTGCGTTGATACCACTGCTTAGATCGGAAGAGCACACGTCTGAA
+
JJJJJIIJJJJJJHIHHHGHFFFFFFCEEEEEDBD?DDDDDDBDDDABDDCA

第一行序列名稱
其中第一行的命名方式在v1.4后是 "@EAS139:\136:\FC706VJ:\2:\2104:\15343:\197393 1:\Y:\18:ATCACG"

在v1.4之前@HWUSI-EAS100R:6:73:941:1973#0/1

** 第三行質(zhì)量序列格式**
目前illumina使用的堿基質(zhì)量格式為phred+33, 和Sanger的質(zhì)量基本一致。

不同版本的堿基質(zhì)量Q和堿基錯誤率P的關(guān)系如下

Relationship between Q and p

FastQC質(zhì)量報告

質(zhì)量控制的軟件很多，但是目前主要以fastqc為主。常見的用法：

fastqc seqfile1 seqfile2 .. seqfileN
常用參數(shù)：
-o： 輸出路徑
--extract: 輸出文件是否需要自動解壓 默認是--noextract
-t: 線程， 和電腦配置有關(guān)，每個線程需要250MB的內(nèi)存
-c: 測序中可能會有污染， 比如說混入其他物種
-a: 接頭
-q: 安靜模式

FastQC有兩種方式分析壓縮的fastq文件

zcat SRR3589956_1.fastq.gz | fastqc -t 4 stdin
fastqc SRR3589956_1.fastq.gz 

結(jié)果會得到一個html文件和一個zip壓縮包。

image.png

其中html文件用瀏覽器打開就能直觀看到數(shù)據(jù)

綠色表示通過，紅色表示未通過，黃色表示不太好。一般而言RNA-Seq數(shù)據(jù)在sequence deplication levels 未通過是比較正常的。畢竟一個基因會大量表達，會測到很多遍。
總體看來，測序可接受。

下面這種（從FASTQC官網(wǎng)找到的實例）就要好好好好處理一下了

bad result

具體含義可以看這里： http://jingyan.baidu.com/article/49711c6149e27dfa441b7c34.html

由于有14個結(jié)果，如果一個一個打開過去，一定會麻煩死，最好有一種一勞永逸的方法。
知乎的青山屋主寫了一篇關(guān)于multiQC的教程（https://zhuanlan.zhihu.com/p/27646873， 介紹聚合多個QC結(jié)果進行演示的方法。

利用conda安裝軟件尤其簡單，

conda install multiqc
multiqc --help

使用也很方便，

# 先獲取QC結(jié)果
 ls *gz | while read id; do fastqc -t 4 $id; done
# multiqc
multiqc *fastqc.zip --pdf

會有一個html文件用來了解總體情況

multiQC

除了用multiQC查看多個QC結(jié)果以外，還可以專門寫一個腳本看每個樣本的reads數(shù)量，GC含量，Q20，Q30的比例

Python腳本，保存為fastqc_summary.py，python3 fastqc_summary.py *.zip，會生成summary.txt文件
邏輯：

1. 用python的zipfile模塊打開zip文件，讀取xx_data.txt的數(shù)據(jù)
2. 讀取每一行的數(shù)據(jù)，用正則表達式模塊re，找到目標行
3. 根據(jù)分隔符對每一行進行分割，進行賦值
4. 由于只需要讀取到>>Per base sequence quality pass這一部分，所以設(shè)置一個>>END_MODULE的計數(shù)器，數(shù)量超過2，就停止。

import re
import zipfile

# read the zip file
def zipReader(file):
    qcfile =  zipfile.ZipFile(file)
    data_txt = [file for file in qcfile.namelist() if re.match(".*?_data\.txt", file)][0]
    data = [bytes.decode(line) for line in qcfile.open(data_txt)]
    return data

def fastqc_summary(data):
    module_num = 0
    bases = 0
    Q20 = 0
    Q30 = 0
    for line in data:
        if re.match('Filename', line):
            filename = line.split(sep="\t")[1].strip()
        if re.match('Total Sequence', line):
            read = line.split(sep="\t")[1].strip()
        if re.match('%GC', line):
            GC = line.split(sep="\t")[1].strip()
        if re.match("[^#](.*?\t){6}",line):
            bases = bases + 1
            if float(line.split("\t")[1]) > 30:
                Q20 = Q20 + 1
                Q30 = Q30 + 1
            elif float(line.split("\t")[1]) > 20:
                Q20 = Q20 + 1

        if re.match(">>END", line) :
            module_num = module_num + 1
            if module_num >= 2:
                break
    Q20 = Q20 / bases
    Q30 = Q30 / bases
    summary = [filename, read, GC, str(Q20), str(Q30)]
    return summary

if __name__ == '__main__':
    import sys
    for arg in range(1, len(sys.argv)):
        data = zipReader(sys.argv[arg])
        summary = fastqc_summary(data)
        with open('summary.txt', 'a') as f:
            f.write('\t'.join(summary) + '\n')