用vs2013做网站案例,广告发布包括哪些,关键词全网搜索,nginx 代理 wordpressLogs
增加R包pctax内的一些帮助上游分析的小脚本#xff08;2024.03.03#xff09;增加Mmseqs2用于去冗余#xff0c;基因聚类的速度非常快#xff0c;且随序列量线性增长#xff08;2024.03.12#xff09;更新全文细节#xff08;2024.05.29#xff09;
注意#x…Logs
增加R包pctax内的一些帮助上游分析的小脚本2024.03.03增加Mmseqs2用于去冗余基因聚类的速度非常快且随序列量线性增长2024.03.12更新全文细节2024.05.29
注意 篇幅有限本文一般不介绍具体的软件安装配置数据库下载等也不提供示例数据主要帮助快速了解整个宏基因组的上游分析流程。所以软件安装请参考软件自身介绍运行时遇到具问题也最好先Google或者去github issue里面问问。 宏基因组分析技术还是在不停的更新迭代很多更新更好用的软件在不断推出所以我在这也会持续更新如果我还一直做这个方向的话。这里介绍的也只是比较基本的分析流程但是掌握了之后自己进一步去做个性化分析也会顺手很多。 完成上游分析后我们可以得到各种物种或者功能的profile后续一般用pythonR进行数据分析和可视化可参考我的其他博文。 绝大多数这里介绍的软件都是仅支持linux平台的我们做测序文件上游分析也肯定是在服务器上做个人PC一般很难满足需求所以在做这些分析前必须先学习linux基础知识如文件系统shell脚本编写软件安装等。安装软件建议使用conda或mamba新建环境和管理有很多参考方法。 我使用的集群是slurm作业管理系统下面的示例脚本也是适用与slurm的尽量先学习一下slurm的使用再尝试提交作业。如果不用slurm的话可以只参考#############中间的软件部分。
Introduction
宏基因组Metagenome是指对一个生态系统中的所有微生物进行DNA分析的过程可以帮助研究人员了解微生物的多样性、功能和互作关系。
宏基因组的应用非常广泛包括
生物多样性研究通过对宏基因组进行分析可以了解不同生态系统中微生物的多样性和分布情况。生态学研究宏基因组可以帮助研究人员了解微生物在生态系统中的功能、互作关系和生态位等。生物技术宏基因组可以用于筛选具有特定功能的微生物例如寻找能够降解有害物质的微生物。
宏基因组的分析一般包括以下步骤
DNA提取与建库。高通量测序使用高通量测序技术对扩增后的DNA进行测序得到原始序列数据。数据清洗和组装对原始数据进行质量控制、去除低质量序列和冗余序列将序列拼接成较长的连续序列contigs。基因注释将contigs中的基因进行注释得到基因功能信息。数据分析了解微生物多样性、群落结构、功能特征等信息更多是指获取了物种丰度表或功能丰度表之后的进一步分析。MAGs binning 进化动态等进一步分析
下图是我常用的一套上游基本流程当然在面对不同项目时应该有不同的侧重点和适用的分析方法可以在此基础上添加或修改。
最早的时候这方面的分析我都是参考刘永鑫老师的EasyMetagenome现在这套流程也发文章了 (1)值得参考对上手16S测序数据或宏基因组数据都很有帮助。 最近2024.3.3我整理了一下流程放在我的R包pctax里了方便整个宏基因组的上下游流程结合 install.packages(pctax)
#使用micro_sbatch可以快速获得每一步的slurm脚本
pctax::micro_sbatch(work_dir ~/work/,step fastp)preprocess
一般把所有样本的测序双端文件放在一个data文件夹下然后把所有样本名写入到一个samplelist文件中方便后续各种软件的批量操作。
质控fastp
测序文件质控是指在下游分析前对测序数据进行质量控制和清理以确保数据的准确性和可靠性。
fastp是一个非常快速、多功能的测序数据质控工具支持自动化的质控处理包括去除低质量序列、剪切接头序列和生成质控报告。
#!/bin/bash
#SBATCH --job-namefastp
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array1
#SBATCH --partitioncpu
#SBATCH --cpus-per-task8
#SBATCH --mem-per-cpu2G
samplelist/share/home/jianglab/pengchen/work/asthma/samplelistecho start: date %Y-%m-%d %T
startdate %s####################
echo SLURM_ARRAY_TASK_ID: $SLURM_ARRAY_TASK_ID
START$SLURM_ARRAY_TASK_IDNUMLINES40 #how many sample in one arraySTOP$((SLURM_ARRAY_TASK_ID*NUMLINES))
START$(($STOP - $(($NUMLINES - 1))))echo START$START
echo STOP$STOP####################
for (( N $START; N $STOP; N ))
dosample$(head -n $N $samplelist | tail -n 1)echo $samplemkdir -p data/fastp~/miniconda3/envs/waste/bin/fastp -w 8 -i data/rawdata/${sample}_f1.fastq.gz -o data/fastp/${sample}_1.fq.gz \-I data/rawdata/${sample}_r2.fastq.gz -O data/fastp/${sample}_2.fq.gz -j data/fastp/${i}.json
done##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $startsfastp的json文件中统计了测序数据的基本指标我们可以将其整理为表格 把所有的.json文件移到一个文件夹里report/下使用我写的R包pctax读取json文件并整理成表格
pctax::pre_fastp(report/)另外FastQCCutadaptTrimmomatic等也是常用的测序质控工具。
去宿主bowtie2
去宿主的过程其实就是将序列比对到宿主基因组上然后没有比对到的序列整合成新文件就是去宿主后的了。宿主基因组需要自己先下载好并用 bowtie2-build 建立索引以人类为例
#!/bin/bash
#SBATCH --job-namerm_human
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array1
#SBATCH --partitioncpu
#SBATCH --cpus-per-task8
#SBATCH --mem-per-cpu2G
samplelist/share/home/jianglab/pengchen/work/asthma/samplelistecho start: date %Y-%m-%d %T
startdate %s####################
echo SLURM_ARRAY_TASK_ID: $SLURM_ARRAY_TASK_ID
START$SLURM_ARRAY_TASK_IDNUMLINES40 #how many sample in one arraySTOP$((SLURM_ARRAY_TASK_ID*NUMLINES))
START$(($STOP - $(($NUMLINES - 1))))echo START$START
echo STOP$STOP####################
for (( N $START; N $STOP; N ))
dosample$(head -n $N $samplelist | tail -n 1)echo $samplemkdir -p data/rm_human~/miniconda3/envs/waste/bin/bowtie2 -p 8 -x ~/db/humangenome/hg38 -1 data/fastp/${sample}_1.fq.gz \-2 data/fastp/${sample}_2.fq.gz -S data/rm_human/${sample}.sam \--un-conc data/rm_human/${sample}.fq --very-sensitiverm data/rm_human/${sample}.sam
done##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $starts另外也可以使用bwakneaddata等软件去宿主。
fastq信息统计
可以用一些小工具如FastqCountseqkit等来统计reads数碱基数GCQ20Q30等指标
~/biosoft/FastqCount-master/FastqCount_v0.5 xx.fastq.gzTotal Reads Total Bases N Bases Q20 Q30 GC
11568822 (11.57 M) 1702829127 (1.70 G) 0.00% 98.00% 94.00% 54.00%还可以把cleandata再跑一下fastp但是只看处理前的统计指标。因为fastp确实非常快纯粹用来统计也行。
reads-based
reads-based宏基因组分析是指直接利用测序读长reads进行宏基因组数据分析的方法而不是先进行基因组组装。该方法通过对短读长进行比对和注释从中提取功能信息和物种组成常用于快速评估环境样本中的微生物群落结构和功能潜力。
优点 快速无需组装处理速度较快。全面能够检测到低丰度的微生物和功能基因适合大规模样本分析。 局限性 片段化由于不进行组装分析基于短读长可能会错过长距离的基因联系和结构信息。依赖数据库分析结果依赖于参考数据库的全面性和准确性。
物种注释kraken2
Kraken2是一个用于对高通量测序数据进行分类和标识物种的软件。它使用参考数据库中的基因组序列来进行分类并使用k-mer方法来实现快速和准确的分类。
使用Kraken2进行基本分类的简单步骤
准备参考数据库Kraken2需要一个参考数据库以便对测序数据进行分类。可以直接下载官方构建的标准库也可以从NCBI、Ensembl或其他数据库下载相应的基因组序列并使用Kraken2内置的工具来构建数据库。
kraken2-build --standard --threads 24 --db ./
--standard标准模式下只下载5种数据库古菌archaea、细菌bacteria、人类human、载体UniVec_Core、病毒viral。 安装Kraken2可以从Kraken2官方网站下载并安装Kraken2软件。 运行Kraken2使用Kraken2对测序数据进行分类需要使用以下命令
kraken2 --db path_to_database input_file --output output_file
这里path_to_database是参考数据库的路径input_file是需要进行分类的输入文件output_file是输出文件的名称。Kraken2将输出一个分类报告文件和一个序列文件。
需要注意的是kraken运行至少要提供数据库大小的内存大小运行内存因为它会把整个数据库载入内存后进行序列的注释所以如果发现无法载入数据库的报错可以尝试调大内存资源。
#!/bin/bash
#SBATCH --job-namekraken2
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --partitioncpu
#SBATCH --cpus-per-task8
#SBATCH --mem-per-cpu40G
samplelist/share/home/jianglab/pengchen/work/asthma/samplelistecho start: date %Y-%m-%d %T
startdate %s####################
echo SLURM_ARRAY_TASK_ID: $SLURM_ARRAY_TASK_ID
START$SLURM_ARRAY_TASK_IDNUMLINES40 #how many sample in one arraySTOP$((SLURM_ARRAY_TASK_ID*NUMLINES))
START$(($STOP - $(($NUMLINES - 1))))####################
for (( N $START; N $STOP; N ))
dosample$(head -n $N $samplelist | tail -n 1)echo $samplemkdir -p result/kraken2~/miniconda3/envs/waste/bin/kraken2 --threads 8 \--db ~/db/kraken2/stand_krakenDB \--confidence 0.05 \--output result/kraken2/${sample}.output \--report result/kraken2/${sample}.kreport \--paired \data/rm_human/${sample}.1.fq \data/rm_human/${sample}.2.fq
done##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $starts另外kraken数据库是可以自己构建的所以适用于各种项目的物种注释我做的比较多的是环境样本的宏基因组就可能需要更全面的物种数据库甚至除了各种微生物还要动植物数据等我们实验室的WX师姐收集构建了一个超大的物种库。
kraken软件运行时载入数据库是一个十分耗时的步骤而每条序列的鉴定时间差不多所以我们可以将很多样本的fastq文件合并成一个大文件后输入kraken注释之后再按照序列的数量拆分结果文件这样多个样本也只需要载入一次数据库节省时间以下是用自定义的kraken2M.py脚本实现。
#!/bin/bash
#SBATCH --job-namekraken2M
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/kraken/%x_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/kraken/%x_%a.err
#SBATCH --time14-00:00:00
#SBATCH --partitionmem
#SBATCH --cpus-per-task32
#SBATCH --mem-per-cpu100G
echo start: date %Y-%m-%d %T
startdate %s####################mkdir -p result/krakenpython /share/home/jianglab/shared/krakenDB/K2ols/kraken2M.py -t 16 \-i data/rm_human \-c 0.05 \-s .1.fq,.2.fq \-o result/kraken \-d /share/home/jianglab/shared/krakenDB/mydb2 \-k ~/miniconda3/envs/waste/bin/kraken2 \-kt /share/home/jianglab/shared/krakenDB/K2ols/KrakenToolsmkdir -p result/kraken/kreportmv result/kraken/*_kreport.txt result/kraken/kreport/python ~/db/script/format_kreports.py -i result/kraken/kreport \-kt /share/home/jianglab/shared/krakenDB/K2ols/KrakenTools --save-name2taxid##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $startsKraken输出格式
标准输出格式output文件(五列表):
C/U代表分类classified或非分类unclassifed序列ID物种注释比序列注释的区域如98|94代表左端98bp右端94bp比对至数据库LCA比对结果如”562:13 561:4”代表13 k-mer比对至物种#5624 k-mer比对至#561物种
报告输出格式report文件(包括6列方便整理下游分析):
百分比countcount最优(U)nclassified, ®oot, (D)omain, (K)ingdom, §hylum, ©lass, (O)rder, (F)amily, (G)enus, or (S)pecies. “G2”代表位于属一种间NCBI物种ID科学物种名
常用的物种丰度表格式除了kraken report还有mpaspfkrona等格式关于kraken结果的整理以及格式转换方式有一些现成的脚本或者自己写。
KrakenTools (jhu.edu) 就是一套很好用的kraken工具包其中常用的有 extract_kraken_reads.py 此程序提取读取在任何用户指定的分类id处分类的内容。用户必须指定Kraken输出文件、序列文件和至少一个分类法ID。下面指定了其他选项。截至2021年4月19日此脚本与KrakenUniq/Kraken2Uniq报告兼容。 combine_kreports.py 这个脚本将多个Kraken报告合并到一个合并的报告文件中。 kreport2krona.py 这个程序需要一个Kraken报告文件并打印出一个krona兼容的文本文件换成krona文件好画图。
转换后的krona使用 ktImportText MYSAMPLE.krona -o MYSAMPLE.krona.html得到网页可视化结果。 kreport2mpa.py 这个程序需要一个Kraken报告文件并打印出一个mpa (MetaPhlAn)风格的文本文件。 combine_mpa.py 这个程序合并来自kreport2mpa.py的多个输出。
物种功能HUMAnN
HUMAnNThe HMP Unified Metabolic Analysis Network是一款用于分析人类微生物组的功能和代谢能力的工具。
它通过将宏基因组序列与参考基因组数据库比对利用MetaCyc代谢通路数据库和UniRef蛋白质序列数据库分析微生物组在功能和代谢通路水平上的组成和活性。HUMAnN2还提供了多样性分析、关联分析和可视化工具可用于深入研究人类微生物组对宿主健康的影响和治疗策略的制定等方面。
HUMAnN是由美国国家人类微生物组计划HMP开发的目前最新版本为HUMAnN3于2020年发布。与HUMAnN2相比HUMAnN3改进了基因家族注释的方法提高了注释精度和速度并提供了新的功能和工具如功能韧度分析、代谢指纹识别和多样性分析等。
但是HUMAnN的数据库基本都是与人相关的微生物比较适合做各种人体微生物组肠道肺部口腔皮肤等等对于环境样本可能unclassified比较多。
HUMAnN要求双端序列合并的文件作为输入for循环根据实验设计样本名批量双端序列合并。
物种组成调用MetaPhlAn3, bowtie2比对至核酸序列库解决有哪些微生物存在的问题功能组成为humann3调用diamond比对至蛋白库解决这些微生物参与哪些功能通路的问题
cd data/rm_human
for i in cat ~/work/asthma/samplelist
do
echo $i
cat ${i}_f1.fastq ${i}_r2.fastq ${i}_paired.fastq
done#!/bin/bash
#SBATCH --job-namehumann
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array1
#SBATCH --partitioncpu
#SBATCH --nodes1
#SBATCH --tasks-per-node1
#SBATCH --cpus-per-task1
#SBATCH --mem-per-cpu2G
samplelist/share/home/jianglab/pengchen/work/asthma/samplelistecho start: date %Y-%m-%d %T
startdate %s####################
echo SLURM_ARRAY_TASK_ID: $SLURM_ARRAY_TASK_ID
START$SLURM_ARRAY_TASK_IDNUMLINES40 #how many sample in one arraySTOP$((SLURM_ARRAY_TASK_ID*NUMLINES))
START$(($STOP - $(($NUMLINES - 1))))#set the min and max
if [ $START -lt 1 ]
thenSTART1
fi
if [ $STOP -gt 40 ]
thenSTOP40
fiecho START$START
echo STOP$STOP####################
for (( N $START; N $STOP; N ))
dosample$(head -n $N $samplelist | tail -n 1)echo $samplemkdir -p data/pairedcat data/rm_human/${sample}.1.fq data/rm_human/${sample}.2.fq data/paired/${sample}_paired.fqmkdir -p result/humann~/miniconda3/envs/waste/bin/humann --input data/paired/${sample}_paired.fq \--output result/humann/ --threads 24ln result/humann/${sample}_paired_humann_temp/${sample}_paired_metaphlan_bugs_list.tsv result/humann/rm -rf result/humann/${sample}_paired_humann_temp
done##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $startsmkdir -p result/metaphlan2
## 合并、修正样本名
merge_metaphlan_tables2.py \result/humann/*_metaphlan_bugs_list.tsv | \sed s/_metaphlan_bugs_list//g \ result/metaphlan2/taxonomy.tsvcontigs-based
contigs-based宏基因组分析是指先将测序读长reads组装成较长的序列片段contigs然后对这些contigs进行进一步的分析。这种方法可以提供更长的基因组片段有助于更准确地进行基因注释和微生物群落结构分析。
优点 更高的解析度通过组装生成较长的序列片段可以更准确地进行基因和基因组注释。结构信息能够获得长距离的基因联系和基因组结构信息有助于更全面的功能和进化分析。 局限性 计算资源需求高组装过程需要较高的计算资源和时间。复杂性组装和binning步骤增加了分析的复杂性需要更多的技术经验和工具支持。
组装assembly
组装assembly是指将测序产生的短读长reads拼接成更长的连续序列contigs的过程。这个过程在基因组学和宏基因组学研究中至关重要因为它能够将片段化的信息整合成更完整的基因组序列便于后续的功能和分类分析。
注意assembly有单样本拼装和混拼等策略可以自行根据计算资源和研究目的选择。
单样本拼装Single-sample Assembly
定义单样本拼装是指对来自单一环境或条件下的样本进行组装。这种策略适用于相对简单的样本或者希望单独分析每个样本的基因组组成。
优点
独立分析每个样本单独组装有助于深入分析样本特定的基因组特征。数据简单适用于数据复杂度较低的样本组装结果较为清晰。
缺点
资源需求高对每个样本进行独立组装计算资源和时间需求较高。有限的覆盖度对于低丰度微生物单样本组装可能无法提供足够的覆盖度和组装质量。
适用场景 环境样本较为简单的研究。 需要对每个样本进行独立比较和分析的研究。 混拼策略Co-assembly
定义混拼是指将多个样本的数据合并在一起进行组装。这种策略适用于复杂的宏基因组样本通过整合多个样本的数据可以提高组装的覆盖度和质量。
优点
提高覆盖度合并多个样本的数据增加了序列的覆盖度有助于更完整的组装。处理复杂样本适用于复杂的环境样本能够捕捉到更多的低丰度微生物和基因组信息。
缺点
复杂性增加混拼的数据复杂度高组装和后续分析更为复杂。样本间差异模糊化样本间的独特特征可能在混拼过程中被模糊化影响对个体样本的具体分析。
适用场景
环境样本复杂包含多种微生物群落。需要提高低丰度物种的检测能力。目标是获取整体微生物群落的综合信息。
megahit
Megahit是一个用于对高通量测序数据进行de novo组装的软件。 它使用了一种基于短读比对和图形构建的算法来组装基因组可以高效地处理大规模的数据集。 以下是Megahit的一些优点和适用情况
速度快Megahit的算法非常高效可以处理大规模的数据集通常比其他de novo组装工具更快。高质量的组装Megahit在组装结果的连通性和准确性方面表现优异尤其在处理高GC含量基因组时效果显著。适用于不同类型的测序数据Megahit支持多种不同类型的测序数据包括 Illumina HiSeq/MiSeq、IonTorrent和PacBio等平台。易于使用Megahit具有简单的命令行语法方便用户进行组装操作且具有中断点避免失败后全部重跑。
#!/bin/bash
#SBATCH --job-namemegahit
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array1
#SBATCH --partitioncpu
#SBATCH --cpus-per-task8
#SBATCH --mem-per-cpu10G
samplelist/share/home/jianglab/pengchen/work/asthma/samplelistecho start: date %Y-%m-%d %T
startdate %s####################
echo SLURM_ARRAY_TASK_ID: $SLURM_ARRAY_TASK_ID
START$SLURM_ARRAY_TASK_IDNUMLINES40 #how many sample in one arraySTOP$((SLURM_ARRAY_TASK_ID*NUMLINES))
START$(($STOP - $(($NUMLINES - 1))))echo START$START
echo STOP$STOP####################
for (( N $START; N $STOP; N ))
dosample$(head -n $N $samplelist | tail -n 1)echo $sample#single samplemkdir -p result/megahitmkdir -p result/megahit/contigs~/miniconda3/envs/waste/bin/megahit -t 8 -1 data/rm_human/${sample}.1.fq \-2 data/rm_human/${sample}.2.fq \-o result/megahit/${sample} --out-prefix ${sample}sed -i //s//${sample}_/ result/megahit/${sample}/${sample}.contigs.famv result/megahit/${sample}/${sample}.contigs.fa result/megahit/contigs/
done##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $starts混拼的话
#!/bin/bash
#SBATCH --job-namemegahit2
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array1
#SBATCH --partitionmem
#SBATCH --cpus-per-task8
#SBATCH --mem-per-cpu100G
samplelist/share/home/jianglab/pengchen/work/asthma/samplelistecho start: date %Y-%m-%d %T
startdate %s#####################multi-sample混拼#需要大内存建议3倍fq.gz内存mkdir -p data/com_readmkdir -p result/megahit2for i in cat samplelistdoecho ${i}cat data/rm_human/${i}.1.fq data/com_read/R1.fqcat data/rm_human/${i}.2.fq data/com_read/R2.fqdone~/miniconda3/envs/waste/bin/megahit -t 8 -1 data/com_read/R1.fq \-2 data/com_read/R2.fq \-o result/megahit2/ --out-prefix multi_sample##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $startsmetaSPAdes
SPAdes 是一款多功能工具包专为测序数据的组装和分析而设计。 SPAdes 主要是为 Illumina 测序数据而开发的但也可用于 IonTorrent。大多数 SPAdes 管道支持混合模式即允许使用长读取PacBio 和 Oxford Nanopore作为补充数据。
而metaSPAdes是目前宏基因组领域组装指标较好的软件尤其在株水平组装优势明显。相比Megahit表现出更好的基因长度和读长比较率但是更好的组装质量也伴随着更长时间和内存消耗同时也有错误组装上升的风险。
#!/bin/bash
#SBATCH --job-nameasthma_megahit
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/megahit/log/%x_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/megahit/log/%x_%a.err
#SBATCH --array1-33
#SBATCH --partitioncpu
#SBATCH --cpus-per-task32echo start: date %Y-%m-%d %T
startdate %s
echo SLURM_ARRAY_TASK_ID: $SLURM_ARRAY_TASK_ID
sample$(head -n $SLURM_ARRAY_TASK_ID ~/work/asthma/data/namelist | tail -1)
#sample$(head -n 1 namelist | tail -1)
echo handling: $sample
####################
for (( N $START; N $STOP; N ))
dosample$(head -n $N $samplelist | tail -n 1)echo $sample#single samplemkdir -p result/metaspadesmkdir -p result/metaspades/contigs~/miniconda3/envs/metawrap/bin/metaspades.py –t 8 -k 21,33,55,77,99,127 --careful \-1 data/rm_human/${sample}.1.fq \-2 data/rm_human/${sample}.2.fq \-o result/metaspades/${sample}sed -i //s//${sample}_/ result/metaspades/${sample}/scaffolds.fastamv result/metaspades/${sample}/scaffolds.fasta result/metaspades/contigs/
done
####################
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $starts组装评估QUAST
QUAST是一个质量评估工具。 QUAST可以使用参考基因组以及不使用参考基因组来评估装配。 QUAST生成详细的报告表格和图解以显示装配的不同方面。
基因预测Prodigal
基因预测是指在基因组序列如contigs中识别出编码区序列CDS, Coding DNA Sequences的过程。这个过程对于理解基因组功能和注释基因组中的基因是至关重要的。常用的软件有ProdigalGeneMarkGeneMarkS用于细菌GeneMark-ES用于真核生物Glimmer等。
这里以Prodigal为例 输入文件拼装好的序列文件 megahit/final.contigs.fa 输出文件prodigal预测的基因序列 prodigal/gene.fa
prodigal不支持多线程运行所以我们可以自行分割序列文件调用多个prodigal程序分别跑实现伪多线程。
#!/bin/bash
#SBATCH --job-nameprodigal
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array1
#SBATCH --partitioncpu
#SBATCH --cpus-per-task1
#SBATCH --mem-per-cpu2G
samplelist/share/home/jianglab/pengchen/work/asthma/samplelistecho start: date %Y-%m-%d %T
startdate %s####################
echo SLURM_ARRAY_TASK_ID: $SLURM_ARRAY_TASK_ID
START$SLURM_ARRAY_TASK_IDNUMLINES40 #how many sample in one arraySTOP$((SLURM_ARRAY_TASK_ID*NUMLINES))
START$(($STOP - $(($NUMLINES - 1))))echo START$START
echo STOP$STOP####################
for (( N $START; N $STOP; N ))
dosample$(head -n $N $samplelist | tail -n 1)echo $samplemkdir -p result/prodigalmkdir -p result/prodigal/gene_famkdir -p result/prodigal/gene_gff~/miniconda3/envs/waste/bin/prodigal -i result/megahit/contigs/${sample}.contigs.fa \-d result/prodigal/gene_fa/${sample}.gene.fa \-o result/prodigal/gene_gff/${sample}.gene.gff \-p meta -f gffmkdir -p result/prodigal/fullgene_fa## 提取完整基因(完整片段获得的基因全为完整如成环的细菌基因组)grep partial00 result/prodigal/gene_fa/${sample}.gene.fa | cut -f1 -d | sed s/// result/prodigal/gene_fa/${sample}.fullidseqkit grep -f result/prodigal/gene_fa/${sample}.fullid result/prodigal/gene_fa/${sample}.gene.fa result/prodigal/fullgene_fa/${sample}.gene.fadone##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $starts去冗余
基因集去冗余是指在预测到的基因集中去除重复的基因或高度相似的基因以获得一个非冗余的基因集合。这一步骤在宏基因组分析中非常重要因为宏基因组样本中通常包含来自多种微生物的重复或相似序列去冗余可以简化数据集并提高后续分析的效率和准确性。常用软件有CD-HITMMseqs2UCLUST集成在QIIME和USEARCH中的工具等。
上面产生了n个样本的基因预测结果文件gene.fa文件要先想办法整合为一个文件再去去冗余。
CD-HIT
之前大家常用的是CD-HIT来去冗余
#!/bin/bash
#SBATCH --job-namecluster
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array1
#SBATCH --partitioncpu
#SBATCH --nodes1
#SBATCH --tasks-per-node1
#SBATCH --cpus-per-task1
#SBATCH --mem-per-cpu2G
samplelist/share/home/jianglab/pengchen/work/asthma/samplelistecho start: date %Y-%m-%d %T
startdate %s####################echo start mergecat result/prodigal/gene_fa/*.gene.fa result/prodigal/all.gene.facat result/prodigal/fullgene_fa/*.gene.fa result/prodigal/all.fullgene.faecho done~/miniconda3/envs/waste/bin/cd-hit-est -i result/prodigal/all.gene.fa \-o result/NR/nucleotide.fa \-aS 0.9 -c 0.9 -G 0 -g 0 -T 0 -M 0mv result/NR/nucleotide_rep_seq.fasta result/NR/nucleotide.farm result/NR/nucleotide_all_seqs.fasta~/miniconda3/envs/waste/bin/seqkit translate result/NR/nucleotide.fa result/NR/protein.fased s/\*//g result/NR/protein.fa result/NR/protein_no_end.farm result/NR/protein.fa##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $startsMmseqs2
最近发现了Mmseqs2这个神器这个聚类要比CD-HIT快非常多其他应用参见Mmseqs2的基础使用
#!/bin/bash
#SBATCH --job-namecluster
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array1
#SBATCH --partitioncpu
#SBATCH --nodes1
#SBATCH --tasks-per-node1
#SBATCH --cpus-per-task1
#SBATCH --mem-per-cpu2G
samplelist/share/home/jianglab/pengchen/work/asthma/samplelistecho start: date %Y-%m-%d %T
startdate %s####################echo start mergecat result/prodigal/gene_fa/*.gene.fa result/prodigal/all.gene.facat result/prodigal/fullgene_fa/*.gene.fa result/prodigal/all.fullgene.faecho donemkdir -p result/NRmmseqs easy-linclust result/prodigal/all.gene.fa result/NR/nucleotide mmseqs_tmp \--min-seq-id 0.9 -c 0.9 --cov-mode 1 --threads 8mv result/NR/nucleotide_rep_seq.fasta result/NR/nucleotide.farm result/NR/nucleotide_all_seqs.fasta~/miniconda3/envs/waste/bin/seqkit translate result/NR/nucleotide.fa result/NR/protein.fased s/\*//g result/NR/protein.fa result/NR/protein_no_end.farm result/NR/protein.fa##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $startsfasta信息统计
使用seqkit统计刚刚得到的一些关键fasta文件的基础信息。
test_filehead -n1 $samplelist
if [ -f result/megahit/contigs/${test_file}.contigs.fa ]; then
~/miniconda3/envs/waste/bin/seqkit stats result/megahit/contigs/*.contigs.fa result/megahit/contig_stats
fi
if [ -f result/prodigal/gene_fa/${test_file}.gene.fa ]; then
~/miniconda3/envs/waste/bin/seqkit stats result/prodigal/gene_fa/*.gene.fa result/prodigal/gene_fa_stats
fi
if [ -f result/prodigal/fullgene_fa/${test_file}.gene.fa ]; then
~/miniconda3/envs/waste/bin/seqkit stats result/prodigal/fullgene_fa/*.gene.fa result/prodigal/fullgene_fa_stats
fi
if [ -f result/NR/nucleotide.fa ]; then
~/miniconda3/envs/waste/bin/seqkit stats result/NR/nucleotide.fa result/NR/nucleotide_stat
fi基因定量salmon
建立索引
#!/bin/bash
#SBATCH --job-namesalmon-index
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array1
#SBATCH --partitioncpu
#SBATCH --cpus-per-task4
#SBATCH --mem-per-cpu10G
echo start: date %Y-%m-%d %T
startdate %s###################### 建索引, -t序列, -i 索引# 大点内存mkdir -p result/salmon~/miniconda3/envs/waste/share/salmon/bin/salmon index \-t result/NR/nucleotide.fa \-p 4 \-i result/salmon/index##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $starts对每个样本定量
#!/bin/bash
#SBATCH --job-namesalmon-quant
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array1
#SBATCH --partitioncpu
#SBATCH --cpus-per-task4
#SBATCH --mem-per-cpu2G
samplelist/share/home/jianglab/pengchen/work/asthma/samplelistecho start: date %Y-%m-%d %T
startdate %s####################
echo SLURM_ARRAY_TASK_ID: $SLURM_ARRAY_TASK_ID
START$SLURM_ARRAY_TASK_IDNUMLINES40 #how many sample in one arraySTOP$((SLURM_ARRAY_TASK_ID*NUMLINES))
START$(($STOP - $(($NUMLINES - 1))))echo START$START
echo STOP$STOP####################
## 输入文件去冗余后的基因和蛋白序列NR/nucleotide.fa
## 输出文件Salmon定量后的结果salmon/gene.count; salmon/gene.TPM
## 定量l文库类型自动选择p线程--meta宏基因组模式
for (( N $START; N $STOP; N ))
dosample$(head -n $N $samplelist | tail -n 1)echo $samplemkdir -p result/salmon/quant~/miniconda3/envs/waste/share/salmon/bin/salmon quant \--validateMappings \-i result/salmon/index -l A -p 4 --meta \-1 data/rm_human/${sample}.1.fq \-2 data/rm_human/${sample}.2.fq \-o result/salmon/quant/${sample}.quant
done##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $starts合并各样本结果
#!/bin/bash
#SBATCH --job-namesalmon-merge
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array1
#SBATCH --partitioncpu
#SBATCH --cpus-per-task1
#SBATCH --mem-per-cpu2G
samplelist/share/home/jianglab/pengchen/work/asthma/samplelistecho start: date %Y-%m-%d %T
startdate %s####################ls result/salmon/quant/*.quant/quant.sf |awk -F/ {print $(NF-1)} |sed s/.quant// tmp_finisheddiff_rows -f samplelist -s tmp_finished tmp_diff# 计算结果的行数line_count$( cat tmp_diff| wc -l)# 检查行数是否大于1如果是则终止脚本if [ $line_count -gt 1 ]; thenecho sf文件和samplelist数量不一致exit 1fi##mapping ratescat result/salmon/quant/*.quant/logs/salmon_quant.log |grep Mapping rate |awk {print $NF} tmppaste samplelist tmp result/salmon/mapping_rates## combine~/miniconda3/envs/waste/share/salmon/bin/salmon quantmerge \--quants result/salmon/quant/*.quant \-o result/salmon/gene.TPM~/miniconda3/envs/waste/share/salmon/bin/salmon quantmerge \--quants result/salmon/quant/*.quant \--column NumReads -o result/salmon/gene.countsed -i 1 s/.quant//g result/salmon/gene.*##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $starts功能基因注释
功能基因注释是对非冗余基因集进行分类和功能描述的过程以便理解这些基因在生物体内可能的角色和作用。
上一步已经有了所有的基因和每个样本所有基因的read count定量结果 我们只需要对上一步的基因序列或蛋白质序列进行不同数据库的注释 然后合并注释结果得到的就是功能丰度表。
很多数据库自带软件都是用diamond比对如果没有专用软件的数据库我们也可以自己用diamond比对。
diamond选择–outfmt 6的输出结果和blastp格式一样。
eggNOG (COG/KEGG/CAZy)
EggNOG数据库收集了COGClusters of Orthologous Groups of proteins直系同源蛋白簇,构成每个COG的蛋白都是被假定为来自于一个祖先蛋白因此是orthologs或者是paralogs。通过把所有完整基因组的编码蛋白一个一个的互相比较确定的。在考虑来自一个给定基因组的蛋白时这种比较将给出每个其他基因组的一个最相似的蛋白因此需要用完整的基因组来定义COG这些基因的每一个都轮番地被考虑。如果在这些蛋白或子集之间一个相互的最佳匹配关系被发现那么那些相互的最佳匹配将形成一个COG。这样一个COG中的成员将与这个COG中的其他成员比起被比较的基因组中的其他蛋白更相像。
EggNOG里面已经包含了GOKEGGCAZy等。
## 下载常用数据库注意设置下载位置
mkdir -p ${db}/eggnog5 cd ${db}/eggnog5
## -y默认同意-f强制下载eggnog.db.gz 7.9G4.9G
download_eggnog_data.py -y -f --data_dir ./## 下载方式2(可选)链接直接下载
wget -c http://eggnog5.embl.de/download/emapperdb-5.0.0/eggnog.db.gz ## 7.9G
wget -c http://eggnog5.embl.de/download/emapperdb-5.0.0/eggnog_proteins.dmnd.gz ## 4.9G
gunzip *.gz#!/bin/bash
#SBATCH --job-nameeggnog
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array1
#SBATCH --partitioncpu
#SBATCH --cpus-per-task8
#SBATCH --mem-per-cpu10G
echo start: date %Y-%m-%d %T
startdate %s
###################### 大点内存, 数据库有50G左右source ~/miniconda3/etc/profile.d/conda.shconda activate func## diamond比对基因至eggNOG 5.0数据库, 1~9h默认diamond 1e-3mkdir -p result/eggnogemapper.py --no_annot --no_file_comments --override \--data_dir ~/db/eggnog5 \-i result/NR/protein_no_end.fa \--cpu 8 -m diamond \-o result/eggnog/protein## 比对结果功能注释, 1hemapper.py \--annotate_hits_table result/eggnog/protein.emapper.seed_orthologs \--data_dir ~/db/eggnog5 \--cpu 8 --no_file_comments --override \-o result/eggnog/output## 添表头, 1列为ID9列KO16列CAZy21列COG22列描述sed 1 i Name\tortholog\tevalue\tscore\ttaxonomic\tprotein\tGO\tEC\tKO\tPathway\tModule\tReaction\trclass\tBRITE\tTC\tCAZy\tBiGG\ttax_scope\tOG\tbestOG\tCOG\tdescription \result/eggnog/output.emapper.annotations \ result/eggnog/eggnog_anno_output##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $startsdbCAN2碳水化合物
## dbCAN2 http://bcb.unl.edu/dbCAN2
## 创建数据库存放目录并进入
mkdir -p ${db}/dbCAN2 cd ${db}/dbCAN2
## 下载序列和描述
wget -c http://bcb.unl.edu/dbCAN2/download/CAZyDB.07312020.fa
wget -c http://bcb.unl.edu/dbCAN2/download/Databases/CAZyDB.07302020.fam-activities.txt
## 备用数据库下载地址并解压
#wget -c http://210.75.224.110/db/dbcan2/CAZyDB.07312020.fa.gz
#gunzip CAZyDB.07312020.fa.gz
## diamond建索引
time diamond makedb \--in CAZyDB.07312020.fa \--db CAZyDB.07312020#!/bin/bash
#SBATCH --job-namecazy
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array1
#SBATCH --partitioncpu
#SBATCH --cpus-per-task8
#SBATCH --mem-per-cpu2G
echo start: date %Y-%m-%d %T
startdate %s####################mkdir -p result/dbcan2diamond blastp \--db ~/db/dbcan2/CAZyDB.07312020 \--query result/NR/protein_no_end.fa \--threads 8 -e 1e-5 --outfmt 6 \--max-target-seqs 1 --quiet \--out result/dbcan2/gene_diamond.f6##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $startsCARD抗生素抗性基因ARGs
RGI软件Github主页: https://github.com/arpcard/rgi可以参考之前写的Antibiotics resistance gene。
#!/bin/bash
#SBATCH --job-namergi
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array1
#SBATCH --partitioncpu
#SBATCH --cpus-per-task8
#SBATCH --mem-per-cpu2Gecho start: date %Y-%m-%d %T
startdate %s####################source ~/miniconda3/etc/profile.d/conda.shconda activate rgimkdir -p result/card~/miniconda3/envs/rgi/bin/rgi main \--input_sequence result/NR/protein_no_end.fa \--output_file result/card/protein \--input_type protein --num_threads 8 \--clean --alignment_tool DIAMOND # --low_quality #partial genes##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $startsVFDB毒力因子
官网地址http://www.mgc.ac.cn/VFs/
在官网下载数据库时带有setA 的库为VFDB数据库核心库(set A)而setB为全库(setB), 其中setA仅包含经实验验证过的毒力基因而setB则在setA的基础上增加了预测的毒力基因选择好数据库后直接用blast/diamond即可完成注释。
#!/bin/bash
#SBATCH --job-namevfdb
#SBATCH --output/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.out
#SBATCH --error/share/home/jianglab/pengchen/work/asthma/log/%x_%A_%a.err
#SBATCH --array1
#SBATCH --partitioncpu
#SBATCH --cpus-per-task8
#SBATCH --mem-per-cpu2G
echo start: date %Y-%m-%d %T
startdate %s####################mkdir -p result/vfdbdiamond blastp \--db ~/db/VFDB/VFDB_setB_pro \--query result/NR/protein_no_end.fa \--threads 8 -e 1e-5 --outfmt 6 \--max-target-seqs 1 --quiet \--out result/vfdb/gene_diamond.f6##############
echo end: date %Y-%m-%d %T
enddate %s
echo TIME:expr $end - $starts其他数据库
NCBI NRNon-redundant:包含广泛的已知蛋白序列是常用的参考数据库。Pfam包含蛋白质家族和功能域的信息通过Hidden Markov Models (HMMs) 描述蛋白质功能域。元素循环数据库专注于环境中的元素如碳、氮、硫、磷等的生物地球化学循环相关基因和途径注释。如NCycDBSCycDBPCycDBMCycDB等
这些数据库我们可以自行下载蛋白数据库然后用blast建库比对。
功能注释合并
写一个python脚本将表1基因-功能的对应表与表2基因丰度表合并即不同基因可能注释到相同功能把它们的丰度加在一起得到新表3功能丰度表。
可以参考https://github.com/YongxinLiu/EasyMicrobiome/blob/master/script/summarizeAbundance.py。
mkdir -p result/summ_table
if [ -f result/eggnog/eggnog_anno_output ]; then
# 汇总9列KO16列CAZy按逗号分隔21列COG按字母分隔原始值累加
~/db/script/summarizeAbundance.py \\
-i result/salmon/gene.count \\
-m result/eggnog/eggnog_anno_output \\
-c 9,16,21 -s ,,* -n raw \\
-o result/summ_table/eggnog
sed -i s/^ko:// summ_table/eggnog.KO.raw.txt
fi分箱binning
宏基因组binning是指将不同的序列集合如metagenome序列集合根据它们的物种归类到不同的bins中以便进一步研究它们的组成和功能。这个过程可以将类似的序列组合在一起形成代表不同物种或基因组的bins以便进行后续分析如物种注释、基因组组装等。
以下是常用的宏基因组binning方法 基于聚类的方法该方法使用序列聚类将相似序列分到同一个bin中。一般来说聚类算法可分为两类无监督聚类如k-means、DBSCAN等和有监督聚类如CAMI、MyCC等。 基于组装的方法该方法使用de novo组装来将相似序列组装成连续的序列再根据这些序列的基因组信息来将其分类到不同的bins中。这种方法的优点是可以更好地处理重复序列缺点是需要大量的计算资源和时间。 基于分类器的方法该方法使用机器学习分类器来将序列分配到不同的bins中。这种方法的优点是可以自动学习特征并在处理大规模数据时效率高缺点是需要先建立一个分类器并进行训练。
在进行宏基因组binning时通常需要使用多个方法进行比较以选择最适合数据集的方法。可以使用一些流行的工具来进行binning如MetaBAT、MaxBin、CONCOCT和MEGAN等。这些工具通常包含各种binning方法可以根据数据集和分析目的选择适合的方法。
篇幅限制具体的方法放在另一篇里面讲解吧查看宏基因组分箱流程。
Reference
Y.-X. Liu, Y. Qin, T. Chen, M. Lu, et al., A practical guide to amplicon and metagenomic analysis of microbiome data. Protein Cell. 12, 315–330 (2021). 关注公众号 ‘bio llbug’,获取最新推送。
