fastq&fasta 处理
FASTQ
FASTQ
完整性校验
md5sum
保证文件在传输过程中保持完整，没有丢失内容，一般采用md5校验方式，目前测序公司给定的测序数据都带有md5文件，使用md5sum -c命令检测文件，如果返回OK，说明文件完整
md5sum -c SRR8651554_1.fastq.md5
md5sum -c SRR8651554_5.fastq.md5
文件统计
统计序列条数，碱基总数，reads读长分布等，可以使用seqkit
统计每条序列ATCG四种碱基组成以及质量值分布，可以使用seqtk comp
统计每个位点所有序列ATCG以及质量值分布，seqtk fqchk命令。结果可以绘制碱基质量以及含量分布图。
合并
可以使用seqtk mergerpe进行合并多个fastq文件，其实cat或者zcat也可以合并，不过seqtk的合并方式有一些差别，cat是将一个文件追加到另一个文件结尾，seqtk mergerpe是每次取文件一个单位合并。
过滤短序列
Ion Torrent，pacbio，nanopore测序的fastq文件序列长度并不相同，通常需要过滤较短的序列，例如过滤掉长度小于150bp的序列。可以使用seqtk seq或者seqkit seq进行操作。
转换为列表
可以使用seqkit fx2tb
转换为列表方便根据ID进行处理。
将四行数据转换为一行三列，可以使用awk列表处理程序来进行处理
使用tab2fx将处理好的列表转为fastq格式。
质量值转换
目前测序得到的fastq文件，都采用phred+33的格式，但是如果处理之前的文件，还有可能遇见phred+64的模式，一般软件中包含--phred33或者--phred64选项，当然也可以直接在两种质量值之间进行转换。
QC
fastqc绘制碱基含量分布图与碱基质量分布图，通过这两个图来判断fastq文件质量好坏。如果样品太多可以使用multiqc合并多个结果。
过滤
去adapter
接头adapter主要是指illumina测序时加入的P7接头与P5接头，理论上来说测序从测序引物开始，是测序不到接头的，但是由于部分打断片段或者由于adapter空载，会导致adpter自身连接，以上两种情况都会导致测序reads中包含adapter序列。adapter序列非基因组本身的序列，会干扰分析，因此需要去除掉。一种是直接给定adapter序列，与reads比对，比对上的就把整条reads去掉，另外一种是测序完成之后，给定一个adater list文件，其中包含了含有adapter序列的reads ID列表，给定一个阈值将这些reads去除掉。cutadapt可以根据给定的adapter序列进行过滤。
去除低质量
低质量主要是指去除测序质量错误率较高的位点，一般以Q20作为标准，如果一个碱基质量值低于Q20，则认为一个低质量，如果一条序列中低质量碱基达到一定比例，例如达到40%以上，则过滤掉此条序列。seqtk工具seq功能通过-q，-n可以将低质量碱基进行标记，例如替换为小写字母或者其他字符，但是不进行过滤，有专门的数据过滤工具。
去除N碱基
如果测序仪无法准确判断出测序碱基的类型，则选择输出N，N碱基无法进行各种分析，因此需要去除掉序列中包含过多N碱基的数据。去除N碱基并不是讲N碱基切除，而且去除包含N碱基过多的整条数据，例如N碱基含量达到10%以上，则过滤掉数据，有些程序按照连续N碱基比率进行过滤。
去除Duplication
duplication是指一对完全一样的测序数据，是由于打断不随机或者PCR过程中导致的，duplication会干扰序列拼接，还会影响变异检查，因此要去除掉。但是RNAseq和宏基因组测序由于本身序列短，并且丰度不同，因此不能去除duplication。去除dupilication 数据可以只保留一对数据，去除多余一致的测序片段，但是在变异检测过程中采取的是在bam文件中对比对到同一位置的duplication片段进行标记的方法，称为Mark Duplication。因为比较reads是否为duplication比较消耗资源，而采用标记的方法更加快速。一般duplication与其他过滤条件一起过滤，或者采用比对之后标记的方式。fastx-toolkit工具中可以去除duplicantion，但只能处理单端，因此用处不大。
截取头尾
illumina测序一般头部有些波动，尾部质量较差，如果想截取部分区域，可以使用seqtk trim功能
数据过滤
有很多软件可以一次性完成数据过滤工作，包括去除adapter，去除低质量，去除N碱基，去除duplication，截取reads，常用的包括fastp，trimmomatic，SOAPnuke等，SOAPnuke很好用，但是去除adapter需要提供一个adapter reads ID的列表，从网上下载的数据没有这个
fastp利用固定adapter序列，但是不能去除duplication
trimmomatic选项参数太长，而且也不是很好用。
这里推荐使用fastp。
排序
如果想对fastq格式文件进行排序，可以使用seqkit sort功能
抽样
有时候需要从全部文件中抽取一部分进行分析，因为测序出来的数据本身就是随机分布的，因此，即使从头到尾开始取数据，出来的也是随机的。当然还是可以随机抽样的。seqtk和seqkit工具都提供了抽样功能。
拆分数据
有时候需要将fastq文件拆成多份，或者根据固定模式进行拆分
例如测序时同一个lane的数据根据index进行拆分，
16S测序中，同一个文件中序列根据barcode进行拆分等。
seqtk split与split2可以用来拆分文件
转换为fasta
一些软件只支持fasta格式，例如只有fasta格式才能进行blast比对，将fastq转换为fasta有多种方法，awk都可以，这里使用seqtk和seqkit分别演示一下。
提取序列
fastq格式文件需要使用bgzip压缩，一般的fastq都是采用gzip压缩，需要重新处理文件才行。
合并两条序列
双末端测序的reads来自一条片段的两侧，如果文库大小比较小，测序读长比较长，例如pairend 300，insertsize 500，就有一些片段中间具有overlap区域，因此可以将两条reads进行拼接，连成更长的片段。有利于进行序列拼接，也具有更好的唯一性，例如16S序列分析中常用此步骤。有很多工具例如cope，flash，fastq-join等。
kmer计数
kmer是一段序列，一般将fastq文件按照固定长度（奇数），例如15，从头到尾按照步移数为1开始截取序列，例如1-15,2-16,3-17……然后对kmer频率进行计数，kmer分析可以用来估计基因组大小等，通过kmer频数分布也可以反应测序错误，或者杂合位点的分布。一般认为kmer频数为1的序列包含测序错误位点。可以使用jellyfish对fastq文件进行kmer计数。
fastq到bam
很多分析的第一步就是通过短序列比对将fastq文件转换为bam，包括变异检测，RNAseq等。一些测序仪直接输出bam格式文件，例如Ion Torrent，属于uBam格式。可以将fastq进行短序列比对的工具很多，这里以bwa为例。除了需要fastq格式，还需要fasta格式作为参考序列。
FASTA
image.png
统计
统计序列条数
grep wc
seqkit的统计功能
格式化
调整，例如一行ID一行序列，或者让每行显示固定长度内容。
逐条统计
统计每条序列长度
seqtk grep awk
统计长度并按照长度计算频数
seqtk grep awk sort uniq
成分统计
计算每条序列的长度以及ATCG碱基的组成，seqtk comp
提取序列
根据基因ID，提取序列。
方法一：利用samtools为fasta建立索引，然后提取
方法二：利用seqtk工具，给定一个ID列表
截取序列
如果想截取某条序列固定区域，需要给定一个bed格式ID列表
seqtk subseq
排序
seqkit sort
按照长度过滤
使用seqtk或者seqkit的seq
反向互补
seqtk是一步完成反向互补操作，
seqkit可以单独取反向序列，也可以单独取互补序列。
-r -p
转化大小写
seq -lu
查找重复序列
RepeatMasker
串联重复序列
trf
blast比对
blastn
建立bwa索引
fasta序列可以作为BWA比对的参考序列，比对之前同样创建索引。
翻译成氨基酸
可以使用一些在线工具例如ExPaSy，EMBOSS等。
