数据预处理 Data preprocessing

最后发布时间 : 2023-11-02 23:23:28 浏览量 :

1.1 input文件准备

运行git命令，获得测试数据和样本信息

git clone https://gitee.com/bioinfoFungi/metagenome.git

.
|-- README.md
|-- result
|   `-- metadata.txt
|-- seq
|   |-- C1_1.fq.gz
|   |-- C1_2.fq.gz
|   |-- C2_1.fq.gz
|   `-- C2_2.fq.gz
`-- temp

seq目录中有2个样本Illumina双端测序，4个序列文件；
temp是临时文件夹，存储分析中间文件，结束可全部删除节约空间
result是重要节点文件和整理化的分析结果图表
metadata.txt: 是实验设计和样本信息

head metadata.txt

# SampleID	Group	Replicate	Sex	Individual	GSA	CRR	project_accession
# C1	Cancer	1	Male  	p136	CRA002355	CRR117732	PRJCA002236
# C2	Cancer	2	Male  	p143	CRA002355	CRR117733	PRJCA002236

zcat C1_1.fq.gz | head

# @SRR3586062.883556
# CTTGGGGCTGCTGAGCTTCATGCTCCCCTCCTGCCTCAAGGACAATAAGGAGATCTTCGACAAGCCTGCAGCAGCTCGCATCGACGCCCTCATCGCTGAGG
# +
# CCCFFFFFHHHHHIJJJJJJJIJIJJJJGIJDGIJEIIJIJJJJJJJJIJJJJIJJIJJJJJHHHFFFFECEEEDDDDD?BDDDDDDBDDDDDDDDBBBDD
# @SRR3586062.3376311
# GACGGTGTCCTCAGGACCCTTCAGTGCCTTCATGATCTGCTCAGAGGTGATGGAGTCACGGACGAGATTCGTCGTGTCAGCACGTAGGATGCGGTCGCCTG
# +
# @@@DDDDAFF?DF;EH+ACHIIICHDEHGIGBFE@GCGDGG?D?G@BGHG@FHCGC;CC:;8ABH>BECCBCB>;8ABCCC@A@#################
# @SRR3586062.3139527
# CCTCTATTTCCAACTTGGTCATTTCTCTCTACATATGAAAATGAATTTACACATAAAATAGAAAACATTAATAATAAAATTTTTTTCATTTAACAACTCCT

1.2 FastQC质量评估(可选)

 fastqc seq/*.gz -t 2

1.2 (可选)FastQC质量评估

    # (可选)使用指定位置的(别人安装的)conda
    source /home/liuyongxin/miniconda2/bin/activate
    # 启动软件环境
    conda activate meta
    # 第一次使用软件要记录软件版本，文章方法中必须写清楚
    fastqc --version # 0.11.9
    # time统计运行时间，fastqc质量评估
    # *.gz为原始数据，-t指定多线程
    time fastqc seq/*.gz -t 2

质控报告见seq目录，详细解读请阅读《数据的质量控制软件——FastQC》。

multiqc将fastqc的多个报告生成单个整合报告，方法批量查看和比较

    # 记录软件版本
    multiqc --version # 1.8
    # 整理seq目录下fastqc报告，输出multiqc_report.html至result/qc目录
    multiqc -d seq/ -o result/qc

查看右侧result/qc目录中multiqc_report.html，单击，选择View in Web Browser查看可交互式报告。

1.3 质量控制

    mkdir -p temp/qc

1.3.1 Fastp质量控制环境样品

适用于无宿主污染的环境样品，质控速度快，自动识别接头和低质量，详见：极速的FASTQ文件质控+过滤+校正fastp

    # 单样本质控
    i=C1
    fastp -i seq/${i}_1.fq.gz -o temp/qc/${i}_1.fastq -I seq/${i}_2.fq.gz -O temp/qc/${i}_2.fastq
 
    # 多样本并行
    tail -n+2 result/metadata.txt|cut -f1|rush -j 2 \
      "fastp -i seq/{1}_1.fq.gz -o temp/qc/{1}_1.fastq -I seq/{1}_2.fq.gz -O temp/qc/{1}_2.fastq"

1.3.2 KneadData质控和去宿主

kneaddata是流程，它主要依赖trimmomatic质控和去接头，bowtie2比对宿主，然后筛选非宿主序列用于下游分析。

详细教程和常见问题，阅读：MPB：随机宏基因组测序数据质量控制和去宿主的分析流程和常见问题

    # 记录核心软件版本
    kneaddata --version # 0.7.4
    trimmomatic -version # 0.39
    bowtie2 --version # 2.3.5.1
    # 可只选一行中部分代码点击Run，如选中下行中#号后面命令查看程序帮助
    # kneaddata -h # 显示帮助

检查点：zless/zcat查看可压缩文件，检查序列质量格式(质量值大写字母为标准Phred33格式，小写字母为Phred64，需参考附录：质量值转换)；检查双端序列ID是否重复，如果重名需要在质控前改名更正。参考附录，质控kneaddata，去宿主后双端不匹配——序列改名。

    # 设置某个样本名为变量i，以后再无需修改
    i=C1
    # zless查看压缩文件，空格翻页，按q退出。
    zless seq/${i}_1.fq.gz | head -n4
    # zcat显示压缩文件，head指定显示行数
    zcat seq/${i}_2.fq.gz | head -n4

"|" 为管道符，上一个命令的输出，传递给下一个命令做输入
gzip: stdout: Broken pipe：管道断开。这里是人为断开，不是错误
运行过程中需要仔细阅读屏幕输出的信息

如果序列双端名称一致，改名参见下方代码

(可选) 序列改名，解决NCBI SRA数据双端ID重名问题，详见《MPB：随机宏基因组测序数据质量控制和去宿主的分析流程和常见问题》。

    gunzip seq/*.gz
    sed -i '1~4 s/$/\\1/g' seq/*_1.fq
    sed -i '1~4 s/$/\\2/g' seq/*_2.fq
    # 再次核对样本是否标签有重复
    head seq/C2_1.fq
    head seq/C2_2.fq
    # 结果压缩节省空间
    gzip seq/*.fq
    # pigz是并行版的gzip，没装可使用为gzip
    # pigz seq/*.fq

(可选)单样品质控

若一条代码分割在多行显示时，最好全部选中运行，多行分割的代码行末有一个 “\” 。多行注释命令运行，可全选，按Ctrl+Shift+C进行注释的取消和添加。

以metadata中C1样品本质控为例

输入文件：双端FASTQ测序数据，提供给参数-i，seq/ ${i}_1.fq.gz和 seq/$ _2.fq.gz
参考数据库：宿主基因组索引 -db $/kneaddata/human_genome/hg37dec_v0.1
输出文件：质控后的FASTQ测序数据，在目录temp/qc下面， ${i}_1_kneaddata_paired_1.fastq和$ _1_kneaddata_paired_1.fastq，用于后续分析
软件位置：conda env list查看软件安装位置，请务必根据自己软件和数据库安装位置，修改软件trimmomatic和接头文件位置。

(可选)手动设置trimmomatic程序和接头位置

    程序目录：${soft}/envs/meta/share/trimmomatic/
    # 查看multiqc结果中接头污染最严重的C2_1样本，再到fastqc报告中查看接头序列，复制前20个碱基检索确定接头文件
    grep 'AGATCGGAAGAGCGTCGTGTAGGGAAA' ${soft}/envs/meta/share/trimmomatic/adapters/*
    # 根据实际情况选择单端SE或双端PE，与原序列比较确定为 TruSeq2-PE.fa，目前多为TruSeq3-PE-2.fa，更准确的是问测序公司要接头文件

100万条序列8线程质控和去宿主，耗时~2m。

    i=C1
    # kneaddata位于liuyongxin中的meta环境，基因组名称为Homo_sapiens或hg37dec_v0.1
    # soft=/home/liuyongxin/miniconda2/
    kneaddata -i seq/${i}_1.fq.gz -i seq/${i}_2.fq.gz \
      -o temp/qc -v -t 8 --remove-intermediate-output \
      --trimmomatic ${soft}/envs/meta/share/trimmomatic/ \
      --trimmomatic-options "ILLUMINACLIP:${soft}/envs/meta/share/trimmomatic/adapters/TruSeq2-PE.fa:2:40:15 SLIDINGWINDOW:4:20 MINLEN:50" \
      --reorder --bowtie2-options "--very-sensitive --dovetail" \
      -db ${db}/kneaddata/human_genome/hg37dec_v0.1

    # 查看质控后的结果文件大小，确保不是0
    ls -shtr temp/qc/${i}_1_kneaddata_paired_?.fastq

多样品并行质控

方法1. rush并行管理：注释修改trimmomatic为绝对路径，即修改/home/liuyongxin/miniconda2为你设置的soft变量完整路径，自己用户安装的soft为~/miniconda3(注：存在单引号的代码内不支持变量)

    tail -n+2 result/metadata.txt|cut -f1|rush -j 2 \
      "kneaddata -i seq/{1}_1.fq.gz -i seq/{1}_2.fq.gz \
      -o temp/qc -v -t 3 --remove-intermediate-output \
      --trimmomatic $soft/envs/meta/share/trimmomatic/ \
      --trimmomatic-options 'ILLUMINACLIP:/home/liuyongxin/miniconda2/envs/meta/share/trimmomatic/adapters/TruSeq2-PE.fa:2:40:15 SLIDINGWINDOW:4:20 MINLEN:50' \
      --reorder --bowtie2-options '--very-sensitive --dovetail' \
      -db ${db}/kneaddata/human_genome/hg37dec_v0.1"

(可选)方法2. parallel并行质控和去宿主

    # 记录软件版本
    parallel --version # 20160222
    # 打will cite承诺引用并行软件parallel
    parallel --citation 
    
    # parallel软件说明和使用实例
    # 根据样本列表`:::`并行处理，并行j=2个任务，每个任务t=3个线程，2~7m
    # 运行下面这行，体会下parallel的工作原理
    # ::: 表示传递参数；第一个::: 后面为第一组参数，对应于{1};
    # 第二个::: 后面为第二组参数，对应于{2}，依次替换
    parallel -j 3 --xapply "echo {1} {2}" ::: seq/*_1.fq.gz ::: seq/*_2.fq.gz
    # --xapply保持文件成对，否则将为两两组合，显示如下：
    parallel -j 2 "echo {1} {2}" ::: seq/*_1.fq.gz ::: seq/*_2.fq.gz
    # 从文件列表使用
    parallel -j 3 --xapply "echo seq/{1}_1.fq.gz seq/{1}_2.fq.gz" ::: `tail -n+2 result/metadata.txt|cut -f1`
    
    # 单样本运行成功，且参数设置绝对路径。出现错误`Unrecognized option: -d64`参考**附录，质控Kneaddata，Java版本不匹配——重装Java运行环境**。
    # 每步分析产生多个文件时新建子文件夹
    # 每个线程处理百万序列约10分钟，多线程可加速 j x t 倍
    # 注意此处引物文件必须填写绝对路径，否则无法使用
    time parallel -j 2 --xapply \
      "kneaddata -i seq/{1}_1.fq.gz \
      -i seq/{1}_2.fq.gz \
      -o temp/qc -v -t 3 --remove-intermediate-output \
      --trimmomatic /home/liuyongxin/miniconda2/envs/meta/share/trimmomatic/ \
      --trimmomatic-options 'ILLUMINACLIP:/home/liuyongxin/miniconda2/envs/meta/share/trimmomatic/adapters/TruSeq2-PE.fa:2:40:15 SLIDINGWINDOW:4:20 MINLEN:50' \
      --reorder --bowtie2-options '--very-sensitive --dovetail' \
      -db ${db}/kneaddata/human_genome/Homo_sapiens" \
      ::: `tail -n+2 result/metadata.txt|cut -f1`

质控结果改名、临时文件删除和统计

大文件清理，高宿主含量样本可节约>90%空间

    rm -rf temp/qc/*contam* temp/qc/*unmatched*  temp/qc/*.fq
    ls -l temp/qc/

结果文件链接为新名：awk的system命令批处理系统命令，s为软链(快捷方式)、f为强制(force)

    awk '{system("ln -sf `pwd`/temp/qc/"$1"_1_kneaddata_paired_1.fastq temp/qc/"$1"_1.fastq")}' <(tail -n+2 result/metadata.txt)
    awk '{system("ln -sf `pwd`/temp/qc/"$1"_1_kneaddata_paired_2.fastq temp/qc/"$1"_2.fastq")}' <(tail -n+2 result/metadata.txt)
    ls -l temp/qc/

质控结果汇总

    # 采用kneaddata附属工具kneaddata_read_count_table
    kneaddata_read_count_table --input temp/qc \
      --output temp/kneaddata.txt
    # 筛选重点结果列
    cut -f 1,2,4,12,13 temp/kneaddata.txt | sed 's/_1_kneaddata//' > result/qc/sum.txt
    cat result/qc/sum.txt

    # 用R代码统计下质控结果，可在本地运行
    Rscript -e "data=read.table('result/qc/sum.txt', header=T, row.names=1, sep='\t'); summary(data)"
    # R转换宽表格为长表格
    Rscript -e "library(reshape2); data=read.table('result/qc/sum.txt', header=T,row.names=1, sep='\t'); write.table(melt(data), file='result/qc/sum_long.txt',sep='\t', quote=F, col.names=T, row.names=F)"
    cat result/qc/sum_long.txt
    # 可用 http://www.ehbio.com/ImageGP/ 绘图展示

1.4 (可选)质控后质量评估

整理bowtie2, trimmomatic, fastqc报告，接头和PCR污染率一般小于1%。结果见：result/qc/multiqc_report_1.html

    # 1-2m
    fastqc temp/qc/*_1_kneaddata_paired_*.fastq -t 2
    multiqc -d temp/qc/ -o result/qc/
    # v1.7以后开始使用Python3，v1.8+缺少bowtie2比对结果的统计

基于reads的分析