分析流程 Analysis pipeline

最后发布时间 : 2023-06-05 10:09:12 浏览量 :

质控KneadData

双端序列质控后是否配对的检查

双端序列质控后序列数量不一致是肯定出错了。但即使序列数量一致，也可能序列不对。在运行metawrap分箱时会报错。可以kneaddata运行时添加--reorder来尝试解决。以下提供了检查双端序列ID是否配对的比较代码

    # 文件
    i=sample1
    seqkit seq -n -i temp/qc/${i}_1_kneaddata_paired_1.fastq|cut -f 1 -d '/' | head > temp/header_${i}_1
    seqkit seq -n -i temp/qc/${i}_1_kneaddata_paired_2.fastq|cut -f 1 -d '/' | head > temp/header_${i}_2
    cmp temp/header_${i}_?

Perl环境不匹配

报错'perl binaries are mismatched'的解决

    e=~/miniconda3/envs/meta
    PERL5LIB=${e}/lib/5.26.2:${e}/lib/5.26.2/x86_64-linux-thread-multi

Java不匹配——重装Java运行环境

若出现错误 Unrecognized option: -d64，则安装java解决：

    conda install -c cyclus java-jdk

读长分析HUMAnN2

HUMAnN2减少输出文件加速

HUMAnN2是计算非常耗时的步骤，如果上百个10G+的样本，有时需要几周至几月的分析。以下介绍两种快速完成分析，而且结果变化不大的方法。替换下面for循环为原文中的“双端合并为单个文件”部分代码

方法1. 软件分析不考虑双端信息，只用一端可获得相近结果，且速度提高1倍。链接质控结果左端高质量至合并目录

    for i in `tail -n+2 result/metadata.txt|cut -f1`;do 
      ln -sf `pwd`/temp/qc/${i}_1_kneaddata_paired_1.fastq temp/concat/${i}.fq
    done

方法2. 控制标准样比对时间。测序数据量通常为6~50G，同一样本分析时间可达10h~100h，严重浪费时间而浪费硬盘空间。
可用head对单端分析截取20M序列，即3G，则为80M行

    for i in `tail -n+2 result/metadata.txt|cut -f1`;do 
       head -n80000000 temp/qc/${i}_1_kneaddata_paired_1.fastq  > temp/concat/${i}.fq
    done

metaphlan2无法找到数据库

正常在首次运行时，会自动下载数据库。有时会失败，解决方法：

方法1. 使用软件安装的用户运行一下程序即可下载成功

方法2. 将我们预下载好的数据索引文件,链接到软件安装目录

    db=~/db
    soft=~/miniconda2
    mkdir -p ${soft}/bin/db_v20
    ln -s ${db}/metaphlan2/* ${soft}/bin/db_v20/
    mkdir -p ${soft}/bin/databases
    ln -s ${db}/metaphlan2/* ${soft}/bin/databases/

CRITICAL ERROR: Can not call software version for bowtie2

解决问题思路：

查看文件位置是否处在conda环境中：type bowtie2。如果不在需要手动设置环境变量的顺序，如果位置正确如在(~/miniconda2/envs/humann2/bin/bowtie2)，请往下看；

检测bowtie2运行情况：bowtie2 -h，报错wd.c: loadable library and perl binaries are mismatched (got handshake key 0xde00080, needed 0xed00080)。错误原因为Perl库版本错误，检查Perl库位置：echo $PERL5LIB，错误原因没有指向环境，并手动修改perl库位置

    # 设置你环境变量位置，最好用绝对路径
    e=~/miniconda2/envs/humann2
    PERL5LIB=${e}/lib/5.26.2:${e}/lib/5.26.2/x86_64-linux-thread-multi

metaphlan_hclust_heatmap.py报错AttributeError: Unknown property axisbg

在网上搜索，axisbg和axis_bgcolor为过时的函数，新版为facecolor，修改为新名称即可 (参考：https://blog.csdn.net/qq_41185868/article/details/81842971)

    # 定位文件绝对路径
    file=`type metaphlan_hclust_heatmap.py|cut -f 2 -d '('|sed 's/)//'`
    # 替换函数名称为新版
    sed -i 's/axisbg/facecolor/g' $file

metaphlan2-共有或特有物种网络图

    awk 'BEGIN{OFS=FS="\t"}{if(FNR==1) {for(i=9;i<=NF;i++) a[i]=$i; print "Tax\tGroup"} \
       else {for(i=9;i<=NF;i++) if($i>0.05) print "Tax_"FNR, a[i];}}' \
       result/metaphlan2/taxonomy.spf > result/metaphlan2/taxonomy_highabundance.tsv
       
    awk 'BEGIN{OFS=FS="\t"}{if(FNR==1) {print "Tax\tGrpcombine";} else a[$1]=a[$1]==""?$2:a[$1]$2;}END{for(i in a) print i,a[i]}' \
       result/metaphlan2/taxonomy_highabundance.tsv > result/metaphlan2/taxonomy_group.tsv
    
    cut -f 2 result/metaphlan2/taxonomy_group.tsv | tail -n +2 | sort -u >group
    
    for i in `cat group`; do printf "#%02x%02x%02x\n" $((RANDOM%256)) $((RANDOM%256)) $((RANDOM%256)); done >colorcode
    
    paste group colorcode >group_colorcode
    
    awk 'BEGIN{OFS=FS="\t"}ARGIND==1{a[$1]=$2;}ARGIND==2{if(FNR==1) {print $0, "Grpcombinecolor"} else print $0,a[$2]}' \
       group_colorcode result/metaphlan2/taxonomy_group.tsv > result/metaphlan2/taxonomy_group2.tsv
    
    awk 'BEGIN{OFS=FS="\t"}{if(FNR==1) {print "Tax",$1,$2,$3,$4, $5, $6, $7, $8 } else print "Tax_"FNR, $1,$2,$3,$4, $5,$6, $7, $8}' \
       result/metaphlan2/taxonomy.spf > result/metaphlan2/taxonomy_anno.tsv

生物标志鉴定LEfSe

lefse-plot_cladogram.py：Unknown property axis_bgcolor

若出现错误 Unknown property axis_bgcolor，则修改lefse-plot_cladogram.py里的ax_bgcolor替换成facecolor即可。

    # 查看脚本位置，然后使用RStudio或Vim修改
    type lefse-plot_cladogram.py

物种分类Kraken2

合并样本为表格combine_mpa.py

krakentools中combine_mpa.py，需手动安装脚本，且结果还需调整样本名

    combine_mpa.py \
      -i `tail -n+2 result/metadata.txt|cut -f1|sed 's/^/temp\/kraken2\//;s/$/.mpa/'|tr '\n' ' '` \
      -o temp/kraken2/combined_mpa

序列筛选/去宿主extract_kraken_reads.py

提取非植物33090和动物(人)33208序列、选择细菌2和古菌2157

    mkdir -p temp/kraken2_qc
    parallel -j 3 \
      "/db/script/extract_kraken_reads.py \
      -k temp/kraken2/{1}.output \
      -r temp/kraken2/{1}.report \
      -1 temp/qc/{1}_1_kneaddata_paired_1.fastq \
      -2 temp/qc/{1}_1_kneaddata_paired_2.fastq \
      -t 33090 33208 --include-children --exclude \
      --max 20000000 --fastq-output \
      -o temp/kraken2_qc/{1}_1.fq \
      -o2 temp/kraken2_qc/{1}_2.fq" \
      ::: `tail -n+2 result/metadata.txt|cut -f1`

组装Megahit

序列长度筛选

megahit默认>200，可选 > 500 / 1000 bp，并统计前后变化；如此处筛选 > 500 bp，序列从15万变为3.5万条，总长度从7M下降到3M

    mv temp/megahit/final.contigs.fa temp/megahit/raw.contigs.fa
    seqkit seq -m 500 temp/megahit/raw.contigs.fa > temp/megahit/final.contigs.fa
    seqkit stat temp/megahit/raw.contigs.fa
    seqkit stat temp/megahit/final.contigs.fa

数据太大导致程序中断

报错信息：126 - Too many vertices in the unitig graph (8403694648 >= 4294967294), you may increase the kmer size to remove tons

解决方法：需要增加k-mer，如最小k-mer改为29，不行继续增加或将数据分批次组装

添加参数： --k-min 29 --k-max 141 --k-step 20

组装MetaSpdades

二三代混合组装

    # 3G数据，耗时3h
    i=SampleA
    time metaspades.py -t 48 -m 500 \
      -1 seq/${i}_1.fastq -2 seq/${i}L_2.fastq \
      --nanopore seq/${i}.fastq \
      -o temp/metaspades_${i}

二三代混合组装OPERA-MS

结果卡在第9步polishing，可添加--no-polishing参数跳过此步；短序列只支持成对文件，多个文件需要cat合并

二三代混合组装

    perl ../OPERA-MS.pl \
        --short-read1 R1.fastq.gz \
        --short-read2 R2.fastq.gz \
        --long-read long_read.fastq \
        --no-ref-clustering \
        --num-processors 32 \
        --out-dir RESULTS

二代组装+三代优化

    perl ~/soft/OPERA-MS/OPERA-MS.pl \
        --contig-file temp/megahit/final.contigs.fa \
        --short-read1 R1.fastq.gz \
        --short-read2 R2.fastq.gz \
        --long-read long_read.fastq \
        --num-processors 32 \
        --no-ref-clustering \
        --no-strain-clustering \
        --no-polishing \
        --out-dir temp/opera

结果可用quast或seqkit stat统计对二代组装的改进效果

基因序列prodigal

序列拆分并行预测基因

(可选)以上注释大约1小时完成1M个基因的预测。加速可将contigs拆分，并行基因预测后再合并。

    # 拆分contigs，按1M条每个文件
    n=10000
    seqkit split result/megahit/final.contigs.fa -s $n
    # 生成拆分文件序列列表
	ls result/megahit/final.contigs.fa.split/final.contigs.part_*.fa|cut -f 2 -d '_'|cut -f 1 -d '.' \
	  > temp/split.list
	# 9线程并行基因预测，此步只用单线程且读写强度不大
	time parallel -j 9 \
	  "prodigal -i result/megahit/final.contigs.fa.split/final.contigs.part_{}.fa \
	  -d temp/gene{}.fa  \
	  -o temp/gene{}.gff -p meta -f gff \
	  > temp/gene{}.log 2>&1 " \
	  ::: `cat temp/split.list`
	# 合并预测基因和gff注释文件
	cat temp/gene*.fa > temp/prodigal/gene.fa
	cat temp/gene*.gff > temp/prodigal/gene.gff

基因去冗余cd-hit

两批基因合并cd-hit-est-2d

cd-hit-est-2d 两批次构建非冗余基因集

A和B基因集，分别有M和N个非冗余基因，两批数据合并后用cd-hit-est去冗余，计算量是(M + N) X (M + N -1)

cd-hit-est-2d比较，只有M X N的计算量

    # 计算B中特有的基因
    cd-hit-est-2d -i A.fa -i2 B.fa -o B.uni.fa \
        -aS 0.9 -c 0.95 -G 0 -g 0 \
        -T 96 -M 0 -d 0
    # 合并为非冗余基因集
    cat A.fa B.uni.fa > NR.fa

cd-hit合并多批基因salmon索引时提示ID重复

	# [error] In FixFasta, two references with the same name but different sequences: k141_2390219_1. We require that all input records have a unique name up to the first whitespace (or user-provided separator) character.
	# 错误解决
	mv temp/NRgene/gene.fa temp/NRgene/gene.fa.bak
	# 15G,2m,4G
	seqkit rename temp/NRgene/gene.fa.bak -o temp/NRgene/gene.fa

基因定量salmon

找不到库文件liblzma.so.0

报错信息：error while loading shared libraries: liblzma.so.0
问题描述：直接运行salmon报告，显示找不到lib库，
解决方法：可使用程序完整路径解决问题，alias salmon="${soft}/envs/metagenome_env/share/salmon/bin/salmon"

基因功能数据库

综合功能注释KEGG描述整理

脚本位于 /db/script 目录，https://www.kegg.jp/kegg-bin/show_brite?ko00001.keg 下载htext，即为最新输入文件 ko00001.keg

    kegg_ko00001_htext2tsv.pl -i ko00001.keg -o ko00001.tsv

抗生素抗性CARD

	# 使用3.1.0和3.1.2均有警告，修改序列名至纯字母数数字也无效
	# WARNING 2021-07-08 08:58:00,478 : Exception : <class 'KeyError'> -> '5141' -> Model(1692) missing in database. Please generate new database.
	# WARNING 2021-07-08 08:58:00,478 : Exception : <class 'KeyError'> -> '5141' -> Model(1692)
	# WARNING 2021-07-08 08:58:00,479 : tetM ---> hsp.bits: 60.8 <class 'float'> ? <class 'str'>

抗生素抗性ResFam

数据库：http://www.dantaslab.org/resfams

参考文献：http://doi.org/10.1038/ismej.2014.106

    mkdir -p temp/resfam result/resfam
    # 比对至抗生素数据库 1m
    time diamond blastp \
      --db ${db}/resfam/Resfams-proteins \
      --query result/NR/protein.fa \
      --threads 9 --outfmt 6 --sensitive \
      -e 1e-5 --max-target-seqs 1 --quiet \
      --out temp/resfam/gene_diamond.f6
    # 提取基因对应抗性基因列表
    cut -f 1,2 temp/resfam/gene_diamond.f6 | uniq | \
      sed '1 i Name\tResGeneID' > temp/resfam/gene_fam.list
    # 统计注释基因的比例, 488/19182=2.5%
    wc -l temp/resfam/gene_fam.list  result/salmon/gene.count 
    # 按列表累计丰度
    summarizeAbundance.py \
      -i result/salmon/gene.TPM \
      -m temp/resfam/gene_fam.list \
      -c 2 -s ',' -n raw \
      -o result/resfam/TPM
    # 结果中添加FAM注释，spf格式用于stamp分析
    awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$4"\t"$3"\t"$2} NR>FNR{print a[$1],$0}' \
      ${db}/resfam/Resfams-proteins_class.tsv  result/resfam/TPM.ResGeneID.raw.txt \
      > result/resfam/TPM.ResGeneID.raw.spf

细菌基因组物种注释GTDB

菌的文件名不要存在非字母数字的符号，否则运行会报错。

	# ERROR: ['BMN5'] are not present in the input list of genome to process，但并无此菌，可能是名称 中存在"-"或"."，替换为i
	# 修改metadata
	sed 's/-/i/;s/\./i/' result/metadatab.txt > result/metadata.txt
    # 修改文件名
    awk 'BEGIN{OFS=FS="\t"}{system("mv temp/antismash/"$1".fna temp/antismash/"$2".fna")ll }' <(paste result/metadatab.txt result/metadata.txt|tail -n+2)

pipleline 软件和数据库安装 Software and database installation