FASTQ
完整性校驗(yàn)
md5sum
保證文件在傳輸過(guò)程中保持完整，沒(méi)有丟失內(nèi)容，一般采用md5校驗(yàn)方式，目前測(cè)序公司給定的測(cè)序數(shù)據(jù)都帶有md5文件，使用md5sum -c命令檢測(cè)文件，如果返回OK，說(shuō)明文件完整
md5sum -c SRR8651554_1.fastq.md5
md5sum -c SRR8651554_5.fastq.md5
文件統(tǒng)計(jì)
統(tǒng)計(jì)序列條數(shù)，堿基總數(shù)，reads讀長(zhǎng)分布等，可以使用seqkit
統(tǒng)計(jì)每條序列ATCG四種堿基組成以及質(zhì)量值分布，可以使用seqtk comp
統(tǒng)計(jì)每個(gè)位點(diǎn)所有序列ATCG以及質(zhì)量值分布，seqtk fqchk命令。結(jié)果可以繪制堿基質(zhì)量以及含量分布圖。
合并
可以使用seqtk mergerpe進(jìn)行合并多個(gè)fastq文件，其實(shí)cat或者zcat也可以合并，不過(guò)seqtk的合并方式有一些差別，cat是將一個(gè)文件追加到另一個(gè)文件結(jié)尾，seqtk mergerpe是每次取文件一個(gè)單位合并。
過(guò)濾短序列
Ion Torrent，pacbio，nanopore測(cè)序的fastq文件序列長(zhǎng)度并不相同，通常需要過(guò)濾較短的序列，例如過(guò)濾掉長(zhǎng)度小于150bp的序列。可以使用seqtk seq或者seqkit seq進(jìn)行操作。
轉(zhuǎn)換為列表
可以使用seqkit fx2tb
轉(zhuǎn)換為列表方便根據(jù)ID進(jìn)行處理。
將四行數(shù)據(jù)轉(zhuǎn)換為一行三列，可以使用awk列表處理程序來(lái)進(jìn)行處理
使用tab2fx將處理好的列表轉(zhuǎn)為fastq格式。
質(zhì)量值轉(zhuǎn)換
目前測(cè)序得到的fastq文件，都采用phred+33的格式，但是如果處理之前的文件，還有可能遇見(jiàn)phred+64的模式，一般軟件中包含--phred33或者--phred64選項(xiàng)，當(dāng)然也可以直接在兩種質(zhì)量值之間進(jìn)行轉(zhuǎn)換。
QC
fastqc繪制堿基含量分布圖與堿基質(zhì)量分布圖，通過(guò)這兩個(gè)圖來(lái)判斷fastq文件質(zhì)量好壞。如果樣品太多可以使用multiqc合并多個(gè)結(jié)果。
過(guò)濾
去adapter
接頭adapter主要是指illumina測(cè)序時(shí)加入的P7接頭與P5接頭，理論上來(lái)說(shuō)測(cè)序從測(cè)序引物開(kāi)始，是測(cè)序不到接頭的，但是由于部分打斷片段或者由于adapter空載，會(huì)導(dǎo)致adpter自身連接，以上兩種情況都會(huì)導(dǎo)致測(cè)序reads中包含adapter序列。adapter序列非基因組本身的序列，會(huì)干擾分析，因此需要去除掉。一種是直接給定adapter序列，與reads比對(duì)，比對(duì)上的就把整條reads去掉，另外一種是測(cè)序完成之后，給定一個(gè)adater list文件，其中包含了含有adapter序列的reads ID列表，給定一個(gè)閾值將這些reads去除掉。cutadapt可以根據(jù)給定的adapter序列進(jìn)行過(guò)濾。
去除低質(zhì)量
低質(zhì)量主要是指去除測(cè)序質(zhì)量錯(cuò)誤率較高的位點(diǎn)，一般以Q20作為標(biāo)準(zhǔn)，如果一個(gè)堿基質(zhì)量值低于Q20，則認(rèn)為一個(gè)低質(zhì)量，如果一條序列中低質(zhì)量堿基達(dá)到一定比例，例如達(dá)到40%以上，則過(guò)濾掉此條序列。seqtk工具seq功能通過(guò)-q，-n可以將低質(zhì)量堿基進(jìn)行標(biāo)記，例如替換為小寫(xiě)字母或者其他字符，但是不進(jìn)行過(guò)濾，有專(zhuān)門(mén)的數(shù)據(jù)過(guò)濾工具。
去除N堿基
如果測(cè)序儀無(wú)法準(zhǔn)確判斷出測(cè)序堿基的類(lèi)型，則選擇輸出N，N堿基無(wú)法進(jìn)行各種分析，因此需要去除掉序列中包含過(guò)多N堿基的數(shù)據(jù)。去除N堿基并不是講N堿基切除，而且去除包含N堿基過(guò)多的整條數(shù)據(jù)，例如N堿基含量達(dá)到10%以上，則過(guò)濾掉數(shù)據(jù)，有些程序按照連續(xù)N堿基比率進(jìn)行過(guò)濾。
去除Duplication
duplication是指一對(duì)完全一樣的測(cè)序數(shù)據(jù)，是由于打斷不隨機(jī)或者PCR過(guò)程中導(dǎo)致的，duplication會(huì)干擾序列拼接，還會(huì)影響變異檢查，因此要去除掉。但是RNAseq和宏基因組測(cè)序由于本身序列短，并且豐度不同，因此不能去除duplication。去除dupilication 數(shù)據(jù)可以只保留一對(duì)數(shù)據(jù)，去除多余一致的測(cè)序片段，但是在變異檢測(cè)過(guò)程中采取的是在bam文件中對(duì)比對(duì)到同一位置的duplication片段進(jìn)行標(biāo)記的方法，稱(chēng)為Mark Duplication。因?yàn)楸容^reads是否為duplication比較消耗資源，而采用標(biāo)記的方法更加快速。一般duplication與其他過(guò)濾條件一起過(guò)濾，或者采用比對(duì)之后標(biāo)記的方式。fastx-toolkit工具中可以去除duplicantion，但只能處理單端，因此用處不大。
截取頭尾
illumina測(cè)序一般頭部有些波動(dòng)，尾部質(zhì)量較差，如果想截取部分區(qū)域，可以使用seqtk trim功能
數(shù)據(jù)過(guò)濾
有很多軟件可以一次性完成數(shù)據(jù)過(guò)濾工作，包括去除adapter，去除低質(zhì)量，去除N堿基，去除duplication，截取reads，常用的包括fastp，trimmomatic，SOAPnuke等，SOAPnuke很好用，但是去除adapter需要提供一個(gè)adapter reads ID的列表，從網(wǎng)上下載的數(shù)據(jù)沒(méi)有這個(gè)
fastp利用固定adapter序列，但是不能去除duplication
trimmomatic選項(xiàng)參數(shù)太長(zhǎng)，而且也不是很好用。
這里推薦使用fastp。
排序
如果想對(duì)fastq格式文件進(jìn)行排序，可以使用seqkit sort功能
抽樣
有時(shí)候需要從全部文件中抽取一部分進(jìn)行分析，因?yàn)闇y(cè)序出來(lái)的數(shù)據(jù)本身就是隨機(jī)分布的，因此，即使從頭到尾開(kāi)始取數(shù)據(jù)，出來(lái)的也是隨機(jī)的。當(dāng)然還是可以隨機(jī)抽樣的。seqtk和seqkit工具都提供了抽樣功能。
拆分?jǐn)?shù)據(jù)
有時(shí)候需要將fastq文件拆成多份，或者根據(jù)固定模式進(jìn)行拆分
例如測(cè)序時(shí)同一個(gè)lane的數(shù)據(jù)根據(jù)index進(jìn)行拆分，
16S測(cè)序中，同一個(gè)文件中序列根據(jù)barcode進(jìn)行拆分等。
seqtk split與split2可以用來(lái)拆分文件
轉(zhuǎn)換為fasta
一些軟件只支持fasta格式，例如只有fasta格式才能進(jìn)行blast比對(duì)，將fastq轉(zhuǎn)換為fasta有多種方法，awk都可以，這里使用seqtk和seqkit分別演示一下。
提取序列
fastq格式文件需要使用bgzip壓縮，一般的fastq都是采用gzip壓縮，需要重新處理文件才行。
合并兩條序列
雙末端測(cè)序的reads來(lái)自一條片段的兩側(cè)，如果文庫(kù)大小比較小，測(cè)序讀長(zhǎng)比較長(zhǎng)，例如pairend 300，insertsize 500，就有一些片段中間具有overlap區(qū)域，因此可以將兩條reads進(jìn)行拼接，連成更長(zhǎng)的片段。有利于進(jìn)行序列拼接，也具有更好的唯一性，例如16S序列分析中常用此步驟。有很多工具例如cope，flash，fastq-join等。
kmer計(jì)數(shù)
kmer是一段序列，一般將fastq文件按照固定長(zhǎng)度（奇數(shù)），例如15，從頭到尾按照步移數(shù)為1開(kāi)始截取序列，例如1-15,2-16,3-17……然后對(duì)kmer頻率進(jìn)行計(jì)數(shù)，kmer分析可以用來(lái)估計(jì)基因組大小等，通過(guò)kmer頻數(shù)分布也可以反應(yīng)測(cè)序錯(cuò)誤，或者雜合位點(diǎn)的分布。一般認(rèn)為kmer頻數(shù)為1的序列包含測(cè)序錯(cuò)誤位點(diǎn)?？梢允褂胘ellyfish對(duì)fastq文件進(jìn)行kmer計(jì)數(shù)。
fastq到bam
很多分析的第一步就是通過(guò)短序列比對(duì)將fastq文件轉(zhuǎn)換為bam，包括變異檢測(cè)，RNAseq等。一些測(cè)序儀直接輸出bam格式文件，例如Ion Torrent，屬于uBam格式?？梢詫astq進(jìn)行短序列比對(duì)的工具很多，這里以bwa為例。除了需要fastq格式，還需要fasta格式作為參考序列。
FASTA
統(tǒng)計(jì)
統(tǒng)計(jì)序列條數(shù)
grep wc
seqkit的統(tǒng)計(jì)功能
格式化
調(diào)整，例如一行ID一行序列，或者讓每行顯示固定長(zhǎng)度內(nèi)容。
逐條統(tǒng)計(jì)
統(tǒng)計(jì)每條序列長(zhǎng)度
seqtk grep awk
統(tǒng)計(jì)長(zhǎng)度并按照長(zhǎng)度計(jì)算頻數(shù)
seqtk grep awk sort uniq
成分統(tǒng)計(jì)
計(jì)算每條序列的長(zhǎng)度以及ATCG堿基的組成，seqtk comp
提取序列
根據(jù)基因ID，提取序列。
方法一：利用samtools為fasta建立索引，然后提取
方法二：利用seqtk工具，給定一個(gè)ID列表
截取序列
如果想截取某條序列固定區(qū)域，需要給定一個(gè)bed格式ID列表
seqtk subseq
排序
seqkit sort
按照長(zhǎng)度過(guò)濾
使用seqtk或者seqkit的seq
反向互補(bǔ)
seqtk是一步完成反向互補(bǔ)操作，
seqkit可以單獨(dú)取反向序列，也可以單獨(dú)取互補(bǔ)序列。
-r -p
轉(zhuǎn)化大小寫(xiě)
seq -lu
查找重復(fù)序列
RepeatMasker
串聯(lián)重復(fù)序列
trf
blast比對(duì)
blastn
建立bwa索引
fasta序列可以作為BWA比對(duì)的參考序列，比對(duì)之前同樣創(chuàng)建索引。
翻譯成氨基酸
可以使用一些在線(xiàn)工具例如ExPaSy，EMBOSS等。
