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)姻报。可以使用jellyfish對(duì)fastq文件進(jìn)行kmer計(jì)數(shù)间螟。
fastq到bam
很多分析的第一步就是通過(guò)短序列比對(duì)將fastq文件轉(zhuǎn)換為bam吴旋，包括變異檢測(cè)，RNAseq等厢破。一些測(cè)序儀直接輸出bam格式文件荣瑟，例如Ion Torrent，屬于uBam格式摩泪“恃妫可以將fastq進(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)建索引透典。
翻譯成氨基酸
可以使用一些在線工具例如ExPaSy晴楔，EMBOSS等。
