下机数据处理：拼接、过滤和去嵌合

下机数据处理：拼接、过滤和去嵌合

参考链接：https://mp.weixin.qq.com/s/aHCMS2yXsAGtmrE8VkDAbg

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-RIanoUPR-1606740390975)(C:\Users\12759\AppData\Roaming\Typora\typora-user-images\image-20201129161922274.png)]

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-of48x4Wd-1606740390984)(C:\Users\12759\AppData\Roaming\Typora\typora-user-images\image-20201129162004285.png)]

数据包含腹泻D和健康H断奶仔猪，有1、3、7、11个时间点，每个时间点有8个样本（D1有6个，H1无H1.6）。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-2C9zaLfR-1606740390986)(C:\Users\12759\AppData\Roaming\Typora\typora-user-images\image-20201129162410489.png)]

混合双端、V3-V4区域测序，00.RawData已经进行了样本拆分、barcode去除和引物切除。每个样本文件夹里有5个文件，第一个extendedfrags.fastq文件是拼接后的序列，raw_1fq.gz和raw_2.fq.gz是未去barcode和引物的双端序列；最后两文件是去掉引物和barcode后的原始数据。

处理过程：先将双端序列进行合并，即reads拼接，用的是flash软件，得到extendedfrags.fastq文件；然后利用qiime1 的split_libraries_fastq.py软件过滤掉低质量序列，即tags过滤或质控，得到.fna文件；再利用vsearch软件进行嵌合体过滤。

01Reads拼接

首先根据Barcode序列和PCR扩增引物序列从下机数据中拆分出各样本数据，截去Barcode和引物序列后使用FLASH软件对每个样本的reads进行拼接，得到的拼接序列为原始Tags数据（Raw Tags）。

FLASH拼接的流程：

a. PE reads比对，找到overlap；

b. 当overlap大于设定的最小overlap值时，执行下面操作：

1）计算overlap长度；

2. 计算错配的数目和overlap的长度两者的比值作为overlap的错配率；

3. 如果计算所得overlap错配率小于现有最优overlap错配率，则将其存为新的最优overlap；

4. 如果错配率和最优overlap一致，计算overlap中所有错配的平均质量值；如果这个平均质量值高于现有最优overlap，则将其存为新的最优overlap；

5）此外，flash软件考虑到 3’端序列质量存在系统性降低趋势，其会根据片段长度在保证PE reads重叠区长度的基础上在3’端对PE reads进行部分截取，这样有利于保证重叠区碱基的质量，提高拼接率。

基本用法
$flash read1.fq read2.fq -m 10 -m 0.1 -o output
主要参数说明（前两个参数最重要）：   
-m 拼接时overlap区的最小长度阈值，默认10bp；   
-x overlap区允许的最大碱基错配比率（最大碱基错配数目/overlap区长度），默认为0.25,微生物组分析中通常设置为0.1。   
-p 碱基质量值类型，64或者33；   
-r reads长度；   
-f 片段长度，也就是测序的文库大小；   
-s 文库的偏差；   
-o 输出文件前缀；   
-t 设置线程数，默认为1，FLASH软件支持多线程，速度快；
FLASH拼接默认输出6个结果文件：   output.extendeFrags.fastq 为拼接后的扩增片段序列文件（最终我们要用到的文件）；   output.flash.log 为日志文件，详细记录了拼接过程中的参数和拼接统计的数据；   output.hist 为拼接后的reads长度的统计信息文件；   output.histogram 为拼接后的reads长度直方图文件；   output.notCombined_1.fastq 为拼接不上的reads1序列文件；   output.notCombined_2.fastq 为拼接不上的reads2序列文件。

代码实现：

将原始序列全部复制到seq目录并重命名：

#*_1/2.fq.gz 【已去除 baraode 和 primer 后得到的序列】
#列出00.RawData目录下的所有文件：
find -type f -print
#复制.fq.gz文件到seq目录下：
find "./" -name "*.fq.gz" |xargs -i cp {} "./seq"
#删除seq目录下包含raw的文件：
ls  *raw*.fq.gz  | xargs rm -fr
#批量修改文件名
rename 'FDMP20H004794-1a_L1_' '' *.fq.gz

#manifest文件获取：
awk 'NR==1{print "sample-id\tforward-absolute-filepath\treverse-absolute-filepath"} NR>1{print $1"\t$PWD/seq/"$1"_1.fq.gz\t$PWD/seq/"$1"_2.fq.gz"}' mapping.txt > manifest

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-S8x2xuHN-1606740390990)(C:\Users\12759\AppData\Roaming\Typora\typora-user-images\image-20201129174459622.png)]

利用flash软件拼接reads，三个小循环。

#flash的shell循环脚本
 #!/bin/bash
 #D1
 for i in 1 2 3 4 5 6
 do
    flash D1.${i}_1.fq.gz D1.${i}_2.fq.gz -m 10 -x 0.1 -o D1.${i}
 done
 #H1
 for i in 1 2 3 4 5 7 8
 do
    flash H1.${i}_1.fq.gz H1.${i}_2.fq.gz -m 10 -x 0.1 -o H1.${i}
 done
 #D H 3 7 11
 for a in D H
     do
         for i in 3 7 11
           do
               for z in 1 2 3 4 5 6 7 8
                 do
                    flash ${a}${i}.${z}_1.fq.gz ${a}${i}.${z}_2.fq.gz -m 10 --max-mismatch-density 0.1 -o ${a}${i}.${z}
                 done
           done
  done

#此时生成的extendedfrags.fastq文件和其他的文件都混杂在seq目录下，将其复制到flashseq目录下方便后续操作
cp *.extendedfrags.fastq ../flashseq

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-GpwZjkqE-1606740390992)(C:\Users\12759\AppData\Roaming\Typora\typora-user-images\image-20201129175326554.png)]

02Tags过滤

对于上步拼接得到的Raw Tags，需要经过严格的过滤处理得到高质量的Tags数据（Clean Tags）。此处我们使用QIIME1软件对拼接数据进行过滤，过滤掉含N较多或含低质量碱基较多的序列; 含低质量碱基较多的序列。

Qiime的质控主要有两个：

a）Tags 截取：将Raw Tags从连续低质量值（默认质量阈值为<=19）碱基数达到设定长度（默认长度值为 3）的第一个低质量碱基位点截断；

b）Tags 长度过滤：Tags 经过截取后得到的 Tags 数据集，进一步过滤掉其中连续高质量碱基长度小于 Tags 长度 75%的Tags。

qiime1 split_libraries_fastq.py参数介绍：   
-i 输入序列文件；    
-b 输入barcode文件，上步产生；    
-m 实验设计，依赖样品barcode列表将序列重复命名    
-q 测序质量阈值，20为99%准确率，30为99.9%准确   
–min_per_read_length_fraction 最小高质量序列比例，0.75即只少有连序75%的碱基质量高于20    
–max_barcode_errors barcode允许的碱基错配个数，0为不允许，即使修改通常也不允许，不信你试试
split_libraries_fastq.py --help查看帮助
注：目前文献中使用 Mothur 和 Qiime 进行数据质控的均比较多，由于 Mothur 无人维护，已经很久不更新了，而 Qiime 是有专人维护，而且不断有更新的，所以使用的为 Qiime软件。

代码实现：

conda activate qiime1进入虚拟环境，使用split_libraries_fastq.py --help查看帮助，

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-EzHeDaIy-1606740390994)(C:\Users\12759\AppData\Roaming\Typora\typora-user-images\image-20201129175948416.png)]

#质控循环脚本
#!/bin/bash
#D1
for i in 1 2 3 4 5 6
   do
       split_libraries_fastq.py -i D1.${i}.extendedFrags.fastq --sample_ids D1.${i} -o fna/D1.${i} -q 19 --barcode_type 'not-barcoded'
   done

#H1
for i in 1 2 3 4 5 7 8
    do
       split_libraries_fastq.py -i H1.${i}.extendedFrags.fastq --sample_ids H1.${i} -o fna/H1.${i} -q 19 --barcode_type 'not-barcoded'
    done

#D、H 3 7 11
for a in D H
    do
      for  i in 3 7 11
           do
              for j in 1 2 3 4 5 6 7 8
                  do
                    split_libraries_fastq.py -i ${a}${i}.${j}.extendedFrags.fastq --sample_ids ${a}${i}.${j} -o fna/${a}${i}.${j} -q 19 --barcode_type 'not-barcoded'
                  done
           done
     done
#此时得到seqs.fna文件，名称一样，分布在flashseq/fna下的子目录中。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-43n9c3sF-1606740390995)(C:\Users\12759\AppData\Roaming\Typora\typora-user-images\image-20201129181110588.png)]

03去嵌合体

嵌合体序列由来自两条或者多条模板链的序列组成，如图所示，在扩增序列X的过程中，在序列延伸阶段，只产生了部分X序列延伸阶段就结束了，在下一轮的PCR反应中，这部分序列作为序列Y的引物接着延伸，扩增就会形成X和Y的嵌合体序列

通常在PCR过程中，大概有1%的几率会出现嵌合体序列，而以在16S/18S/ITS 扩增子测序的分析中，系统相似度极高，嵌合体可达1%-20%，需要去除嵌合体序列。嵌合体的比例与PCR循环数相关，循环数越高，嵌合体比例越高。事实上嵌合体在正常生物体中是不存在的，所以在16S扩增子测序的分析中，需要去除嵌合体序列。

将质控得到的Tags与数据库（16S推荐使用Silva database）进行比对检测嵌合体序列并最终去除其中的嵌合体序列，得到最终的有效数据（Effective Tags），从而用于下游分析。此处，我们使用uchime_ref进行有参去嵌合体，数据库推荐使用又大又全的SILVA最新版本。同时推荐不要基于参考序列去嵌合，因为亲本缺少丰度信息情况下，容易造成假阴性。而de novo去嵌合时，要求亲本的丰度至少是嵌合体的16倍以上，这样可以较少控制阴性率。

基本用法：
vsearch -uchime_ref fastafile -nonchimeras outputfile -db fastafile
部分参数介绍：  
-uchime_ref   有参去嵌合体  
-nonchimeras filename 输出无嵌合体结果文件  
-db 文件名  使用使用–uchime_ref指定数据库fasta文件（eg：小编使用的SILVA_138_SSURef_NR99_tax_silva.fasta，这个需要自行下载）

代码实现：

#去嵌合循环脚本
#!/bin/bash
for i in *
   do
       vsearch -uchime_ref "$i"/seqs.fna -nonchimeras "$i"/"$i".fna -db "/home/lixiyang/renxingda/result/SILVA_138_SSURef_NR99_tax_silva.fasta"
done
#经此步得到去嵌合的序列，即cleanseq（如图D1.1fna）。
#为便于后续操作，将cleanseq全部复制到新的文件夹cleandata下
mkdir cleandata ; find . -name *D*.fna | xargs -i cp {} cleandata
find . -name *H*.fna | xargs -i cp {} cleandata
#然后用cat命令拼接所有的cleanseq,得到（即图中cleanseq.fna）。
cat *.fna > cleanseq.fna

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-UkCzLxBt-1606740390998)(C:\Users\12759\AppData\Roaming\Typora\typora-user-images\image-20201130190927584.png)]

qiime1进行otu聚类

参考链接： https://www.cnblogs.com/xudongliang/p/7205190.html

qiime 本身不提供聚类的算法，它只是对其他聚otu软件的封装，根据聚类软件的算法，分成了3个方向：

de novo： pick_de_novo_otus.py

closed-reference： pick_closed_reference_otus.py

open-reference OTU： pick_open_reference_otus.py

利用pick_open_reference_otus.py脚本进行聚类, 共有6个步骤，其中前4步为OTU 聚类，后2步为产生OTU table 和 聚类的tree。

代码实现：

pick_open_reference_otus.py -i $PWD/cleanseq.fna -r $PWD/refseqs.fna -o $PWD/ucrss_sortmerna_sumaclust/
#-i 后接输入文件；-r后接比对数据库文件，fasta格式, 默认采用的是 greengene/97_otus.fasta；-o后接输出目录
#然后将得到的otu_table_mc2.biom进行格式转换为tsv方便查看
biom convert -i otu_table_mc2.biom -o otu_table.txt --to-tsv
#完成后会得到以下文件

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-F8MC46pw-1606740390999)(C:\Users\12759\AppData\Roaming\Typora\typora-user-images\image-20201130201837081.png)]