1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

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

# C1_2.fq.gz
@SRR3586062.883556
AAGCCAGGACCCGAGTAGTCTCGGCCAGCATAGTACTCCATACGCCAGCCCTTGGGCGCAGCTTCGAGAAGGGCTTTCTCCTCAGCGATGAGGGCGTCGAT
+
CCCFFFFFHHHHHJJHIJEHIJJHIIIIIIJJJHIJGIJIGIJJJJGGCHFEGCHGHGFBCDCDDBD?ADDDDDDDDDDDDCDDDDD@BDCCCBDD<BD><
@SRR3586062.3376311
TGGCGTGCAGTTCAATACCCTCCTCGACACCGTCGTCGTCGTGCGTGTCATCGAGGGTGGCCCCAGCGAGCGTGCAGGCCTGCAGCCAGGCGACCGCATCT
+
=8+=?=)0<+CDDEAEEA<C>FD>F<EEDADDBD)@DA7BADDD8CDA(7.;@A@@>=;:?==??88;>9>>;9-<>>>>A82<AA###############

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

wd=PATH/TO/WORKINGDIR
cd $wd

#### 这里选择在kneaddata环境中安装质量控制所需软件
conda activate kneaddata

# fastqc 0.12.1 kneaddata环境自带，无需下载

# fastp 0.23.2
conda install fastp

# multiqc 1.27
conda install multiqc

# rush 0.6.0
conda install rush

# 将python升级到3.9及以上版本
conda install python=3.9

1
2
3
4
5

# 使用fastqc进行质量评估，-t 1代表设置线程数为1
fastqc seq/*.gz -t 1

# 使用multiqc进行多样本汇总，指定测序样本文件夹和结果输出文件夹
multiqc -d seq/ -o result/qc

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

tail -n +2 metadata.txt | cut -f1 | rush -j 2 \
	"fastp -i seq/{}_1.fq.gz -I seq/{}_2.fq.gz \
	-j temp/qc/{}_fastp.json -h temp/qc/{}_fastp.html \
	-o temp/qc/{}_1.fastq  -O temp/qc/{}_2.fastq \
	> temp/qc/{}.log 2>&1 "

# 提取temp/qc/${i}.log里的关键信息，生成质控的摘要，即每个fastq文件的原始数据条目数量以及滤过后的数据条目数量
echo -e "SampleID\tRaw\tClean" > temp/fastp
	for i in `tail -n+2 metadata.txt | cut -f1`; do
		echo -e -n "$i\t" >> temp/fastp
		grep 'total reads' temp/qc/${i}.log | uniq | cut -f2 -d ':' | tr '\n' '\t' >> temp/fastp
		echo "" >> temp/fastp
 	done
sed -i 's/ //g;s/\t$//' temp/fastp

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

# 创建并进入kneaddata环境
conda create -c bioconda -n kneaddata kneaddata
conda activate kneaddata

# 从NCBI上下载人类参考基因组CRCh38.(3.11Gb)
wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.40_GRCh38.p14/GCF_000001405.40_GRCh38.p14_genomic.fna.gz

# 使用bowtie2生成索引文件，这个步骤需要消耗数小时的时间。对于人类基因组可以直接下载生成好的索引文件，而对于其它小众物种的基因组，则仍需要自行使用bowtie2生成。
bowtie2-build GCF_000001405.40_GRCh38.p14_genomic.fna /mnt/d/metagenome/databases/kneaddata/huma
n --threads 10

# 指定数据库的存放位置
db=/mnt/d/metagenome/databases

# 使用kneaddata去宿主。这里需要对原始测序文件进行处理：同一个样本需要有相同的前缀名称，在名称最后加上“/1”或“/2”分别表示是前向测序还是反向测序。这里通过sed命令将所有{}_1.fastq的每条记录的名称末尾加上“/1”，同时将所有{}_2.fastq的每条记录的名称末尾加上“/2”，之后将处理好的fastq文件传给kneaddata进行过滤。
tail -n+2 metadata.txt | cut -f1 | rush -j 2 \
	"sed '1~4 s/^@\(.*\)$/@\1\/1/' temp/qc/{}_1.fastq > /tmp/{}_1.fastq; \
	sed '1~4 s/^@\(.*\)$/@\1\/2/' temp/qc/{}_2.fastq > /tmp/{}_2.fastq; \
    kneaddata -i1 /tmp/{}_1.fastq -i2 /tmp/{}_2.fastq \
    -o temp/hr --output-prefix {} \
    --bypass-trim --bypass-trf --reorder \
    --bowtie2-options '--very-sensitive --dovetail' \
    -db ${db}/kneaddata/human\
    --remove-intermediate-output -v -t 4; \
    rm /tmp/{}_1.fastq /tmp/{}_2.fastq"

# 结果改名
rename 's/paired_//' temp/hr/*.fastq

# 质控结果汇总，即读取“{}.log”文件
kneaddata_read_count_table --input temp/hr \
	--output temp/kneaddata.txt
	
# 筛选重点结果列，我们只关心1、2、5、6列的信息
cut -f 1,2,5,6 temp/kneaddata.txt | sed 's/_1_kneaddata//' > result/qc/sum.txt

# 查看内容
cat result/qc/sum.txt

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

# 安装biobakery/humann3环境
conda create --name biobakery3 python=3.7
conda activate biobakery3
conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge
conda config --add channels biobakery
conda install humann -c biobakery
db=/mnt/d/metagenome/databases

# 下载三个数据库。下载完毕后会在配置文件中自动修改数据库的存放位置。
# upgrade your pangenome database: 微生物物种核心基因(约15Gb)
humann_databases --download chocophlan full ${db}/humann --update-config yes
# upgrade your protein database: 功能基因diamond索引(约34Gb)
humann_databases --download uniref uniref90_diamond ${db}/humann --update-config yes
# upgrade your annotations database: 输助比对数据库(约3Gb)
humann_databases --download utility_mapping full ${db}/humann --update-config yes


# 在用humann3分析之前，需要将双端测序的两个文件合并为一个文件
for i in `tail -n+2 metadata.txt | cut -f1`; do 
	cat temp/hr/${i}_?.fastq \
	> temp/concat/${i}.fq; done
	
# 在MetaPhlAn4 官方网站下载最新的数据库文件，解压到数据库。然后修改metaphlan_database.cfg文件，将bowtie2db和mpa_pkl参数的值改为数据库文件的路径。如果不写--offline参数，则命令运行时会自动下载。(约20Gb)
# 该命令同时运行了metaphlan和humann命令，生成的关键文件包括物种相对丰度文件、通路丰度文件以及基因家族丰度文件。
tail -n+2 metadata.txt | cut -f1 | rush -j 2 \
	"humann --input temp/concat/{}.fq\
	--output temp/humann3/ \
	--threads 8\
    --metaphlan-options '-t rel_ab --sample_id {} --offline'"

# 将${}_metaphlan_bugs_list.tsv文件移动位置
for i in $(tail -n+2 metadata.txt | cut -f1); do  
	mv temp/humann3/${}_humann_temp/${}_metaphlan_bugs_list.tsv temp/humann3/
done

# 合并表格。注意，合并后tsv文件的第二行恒定为“#metaphlan4”，以便接下来的metaphlan_to_stamp.pl能够正确识别。
merge_metaphlan_tables.py temp/humann3/*_metaphlan_bugs_list.tsv | \
	sed 's/_metaphlan_bugs_list//g' |\
    tail -n+2 | sed '1 s/clade_name/ID/' |\
    sed '2i #metaphlan4'> result/metaphlan4/taxonomy.tsv
      
# metaphlan4较2增加更多unclassified和重复结果，用sort和uniq去除：在https://github.com/LangilleLab/microbiome_helper/blob/master/metaphlan_to_stamp.pl上下载metaphlan_to_stamp.pl文件，放到~/miniconda3/envs/biobakery3/bin文件夹内。
chmod +777 ~/miniconda3/envs/biobakery3/bin/metaphlan_to_stamp.pl
metaphlan_to_stamp.pl result/metaphlan4/taxonomy.tsv \
	|sort -r |\
    uniq > result/metaphlan4/taxonomy.spf
	
# STAMP不支持unclassified，需要过滤掉再使用。使用该文件可进行stamp分析
grep -v 'unclassified' result/metaphlan4/taxonomy.spf > result/metaphlan4/taxonomy2.spf
	
# 将文件的最后一行移动至第一行，以便STAMP软件可以正确识别。
(tail -n 1 "result/metaphlan4/taxonomy2.spf" && head -n -1 "result/metaphlan4/taxonomy2.spf") \
	> result/metaphlan4/taxonomy3.spf

1
2

# https://github.com/baidefeng/EasyMetagenome/blob/master/result12/metaphlan4/taxonomy.spf
# https://github.com/baidefeng/EasyMetagenome/blob/master/result12/metadata.txt

1
2

# 过滤掉Unclassified结果
grep -v 'unclassified' result/metaphlan4/taxonomy.spf > result/metaphlan4/taxonomy2.spf

1
2
3
4
5
6

humann_join_tables --input temp/humann3 \
      --file_name pathabundance \
      --output result/humann3/pathabundance.tsv

# 将样本名后面的后缀去掉，更加美观
sed -i 's/_Abundance//g' result/humann3/pathabundance.tsv

1
2
3
4

humann_renorm_table \
	--input result/humann3/pathabundance.tsv \
	--units relab \
	--output result/humann3/pathabundance_relab.tsv

1
2
3

humann_split_stratified_table \
	--input result/humann3/pathabundance_relab.tsv \
	--output result/humann3/ 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

#### 【关键步骤】给上一步生成的stratified表的第二行加上metadata信息，使之变成pcl文件（*.pcl是接下来分析和画图的关键文件！）
# 提取样品列表
head -n1 result/humann3/pathabundance_relab_stratified.tsv | sed 's/# Pathway/SampleID/' | tr '\t' '\n' > temp/header
# 对应分组，本示例分组为第3列($3)，根据实际情况修改。
awk 'BEGIN{FS=OFS="\t"}NR==FNR{a[$1]=$3}NR>FNR{print a[$1]}' metadata.txt temp/header | tr '\n' '\t'|sed 's/\t$/\n/' > temp/group
# 合成样本、分组+数据
cat <(head -n1 result/humann3/pathabundance_relab_stratified.tsv) temp/group <(tail -n+2 result/humann3/pathabundance_relab_stratified.tsv) \
      > result/humann3/pathabundance_relab_stratified.pcl   
      
#### 对上一步生成的pcl文件，进行表达差异分析与柱状图的绘制
# 使用humann_associate命令进行差异分析，其中需要指定metadata的名称和类型，以及p值。
pcl=result/humann3/pathabundance_relab_stratified.pcl   
humann_associate --input ${pcl} \
	--focal-metadatum Replicate --focal-type categorical \
	--last-metadatum Replicate --fdr 0.05 \
	--output result/humann3/associate.txt
wc -l result/humann3/associate.txt

# 指定感兴趣的通路名称，使用humann_barplot绘制分层条形图
path=PWY-6277
humann_barplot \
	--input ${pcl} --focal-feature ${path} \
	--focal-metadata Replicate --last-metadata Replicate \
	--output result/humann3/barplot_${path}.pdf --sort sum metadata

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

# 下载别人的*.pcl文件
wget -c http://www.imeta.science/github/EasyMetagenome/result/humann2/hmp_pathabund.pcl
cp hmp_pathabund.pcl result/humann3/hmp_pathabund.pcl

# 进行差异分析。这里指定使用Group进行分组
pcl=result/humann3/hmp_pathabund.pcl
humann_associate --input ${pcl} \
	--focal-metadatum Group --focal-type categorical \
	--last-metadatum Group --fdr 0.05 \
	--output result/humann3/associate.txt
wc -l result/humann3/associate.txt

# 选择感兴趣的通路进行绘图
path=P163-PWY
humann_barplot \
	--input ${pcl} --focal-feature ${path} \
	--focal-metadata Group --last-metadata Group \
	--output result/humann3/barplot_${path}.pdf --sort sum metadata 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

# 转换基因家族为KO(uniref90_ko)，可选eggNOG(uniref90_eggnog)或酶(uniref90_level4ec)
for i in `tail -n+2 metadata.txt|cut -f1`;do
  humann_regroup_table \
    -i temp/humann3/${i}_genefamilies.tsv \
    -g uniref90_ko \
    -o temp/humann3/${i}_ko.tsv
done

# 合并，并修正样本名
humann_join_tables \
  --input temp/humann3/ \
  --file_name ko \
  --output result/humann3/ko.tsv
sed -i '1s/_Abundance-RPKs//g' result/humann3/ko.tsv

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

# 创建GraPhlAn环境
conda create -n Graphlan
conda activate Graphlan

# 下载GraphlAn程序
conda install graphlan

# 下载格式转化程序，可以很方便的将LEfSe、MetaPhlAn4结果转化为GraPhlAn的输入文件
conda install export2graphlan

# 需要额外下载python库：biopython
conda install biopython

# 下载R语言核心以运行R脚本
conda install r-base
conda install r-optparse

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

# 下载示例文件
wget -c https://github.com/baidefeng/EasyMetagenome/raw/refs/heads/master/result12/metaphlan4/taxonomy.tsv
mv taxonomy.tsv result/metaphlan4/taxonomy.tsv

# 【关键步骤！】将taxonomy.tsv转化成*.tree.txt与*.annot.txt
export2graphlan.py --skip_rows 1,2 \
	-i result/metaphlan4/taxonomy.tsv \
	--tree temp/merged_abundance.tree.txt \
	--annotation temp/merged_abundance.annot.txt \
	--most_abundant 50 \
	--abundance_threshold 10 \
	--least_biomarkers 3 \
	--annotations 1,2,3,4,5,6 \
	--external_annotations 7
  
# 生成*.xml文件
graphlan_annotate.py --annot temp/merged_abundance.annot.txt \
  temp/merged_abundance.tree.txt  temp/merged_abundance.xml
  
# 输出pdf，这里会输出三个pdf文件，即主图，以及两个图注，后续可以使用Adobe Illustrator进行拼接
graphlan.py temp/merged_abundance.xml result/metaphlan4/graphlan.pdf \
  --external_legends 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

cd result

# 下载tsv2spf.py，并使用该python脚本将tsv文件转化为可供绘图的*.spf文件。
wget -c https://github.com/joybio/graphlan_plot/raw/refs/heads/main/scripts/tsv2spf.py
chmod a+x tsv2spf.py
python ./tsv2spf.py -i metaphlan4/taxonomy.tsv -o metaphlan4/taxonomy_plot.spf

# 【关键步骤！】下载R脚本。这里R脚本依赖R包：optparse，在上文中已经下载过。可以使用"conda install r-optparse"下载
wget -c https://github.com/joybio/graphlan_plot/raw/refs/heads/main/graphlan_plot.r

# 边缘以热图进行展示
Rscript ./graphlan_plot.r \
	-i metaphlan4/taxonomy_plot.spf\
	-d metaphlan4/metadata.txt\
	-c Group\
	-g all\
	-n 100\
	-t heatmap\
	-o metaphlan4/Heat_Results
	
# 边缘以条形图进行展示
Rscript ./graphlan_plot.r \
	-i metaphlan4/taxonomy_plot.spf\
	-d metaphlan4/metadata.txt\
	-c Group\
	-g all\
	-n 100\
	-t bar\
	-o metaphlan4/Bar_Results

cd ..

1
2
3

conda creat -n Lefse
conda activate Lefse
conda install lefse

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# 提取样本行替换为每个样本一行，修改ID为SampleID
head -n1 result/metaphlan4/taxonomy.tsv|\
	tr '\t' '\n'|\
	sed '1 s/ID/SampleID/' >temp/sampleid
# 提取SampleID对应的分组Group(假设为metadata.txt中第二列$2)，替换换行\n为制表符\t，再把行末制表符\t替换回换行
awk 'BEGIN{OFS=FS="\t"}NR==FNR{a[$1]=$2}NR>FNR{print a[$1]}' result/metaphlan4/metadata.txt temp/sampleid|tr '\n' '\t'|sed 's/\t$/\n/' >groupid
# 合并分组和数据(替换表头)
sed -n '1p' result/metaphlan4/taxonomy.tsv > sampleid
cat groupid sampleid <(tail -n+2 result/metaphlan4/taxonomy.tsv) >result/metaphlan4/lefse.txt
# 格式转换为lefse内部格式
lefse_format_input.py result/metaphlan4/lefse.txt \
	temp/input.in -c 1 -s -1 -u 2 -o 1000000

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

# 运行Lefse分析程序，输出结果为input.res
lefse_run.py temp/input.in temp/input.res

# 绘制物种数，标注差异
lefse_plot_cladogram.py temp/input.res \
	result/metaphlan4/lefse_cladogram.pdf --format pdf

# 绘制所有差异features的柱状图
lefse_plot_res.py temp/input.res \
	result/metaphlan4/lefse_res.pdf --format pdf
	
# 指定一个feature，绘制组间对比的图
lefse_plot_features.py -f one --format pdf \
	--feature_name "k__Bacteria.p__Firmicutes" \
	temp/input.in temp/input.res \
	result/metaphlan4/lefse_Firmicutes.pdf

1
2
3
4
5
6
7
8
9

conda create -c bioconda -n kraken2 kraken2
conda activate kraken2

# 创建数据库文件夹
DBNAME=databases/kraken2
mkdir -p $DBNAME

#直接下载标准数据库（约92.3GB）下载完解压至上述文件夹内
wget -c https://genome-idx.s3.amazonaws.com/kraken/k2_standard_20241228.tar.gz

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

# 该步骤需要将80Gb左右的hash表读进内存，故电脑内存小于该值的情况下必须要加上--memory-mapping参数将文件映射至虚拟内存中。（速度可能会减慢！）
for i in `tail -n+2 metadata.txt | cut -f1`; do
	kraken2 --db ${DBNAME}\
		--paired temp/hr/${i}_?.fastq \
		--threads 8 --use-names --report-zero-counts \
		--report temp/kraken2/${i}.report \
 		--output temp/kraken2/${i}.output\
 		--memory-mapping
done

for i in `tail -n+2 result/metadata.txt | cut -f1`; do
	kreport2mpa.py -r temp/kraken2/${i}.report \
	--display-header -o temp/kraken2/${i}.mpa; done

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

mkdir -p result/kraken2
# 输出结果行数相同，但不一定顺序一致，要重新排序
tail -n+2 metadata.txt | cut -f1 | rush -j 1 \
	'tail -n+2 temp/kraken2/{1}.mpa |\
	LC_ALL=C sort |\
	cut -f 2 |\
	sed "1 s/^/{1}\n/" > temp/kraken2/{1}_count '
	
# 提取第一样本品行名为表行名
header=`tail -n 1 metadata.txt | cut -f 1`
tail -n+2 temp/kraken2/${header}.mpa | LC_ALL=C sort | cut -f 1 | \
  sed "1 s/^/Taxonomy\n/" > temp/kraken2/0header_count
head -n3 temp/kraken2/0header_count

# paste合并样本为表格
ls temp/kraken2/*count
paste temp/kraken2/*count > result/kraken2/tax_count.mpa