代码涉及软件、数据库都是绝对路径
source /home/dengqr/miniconda3/bin/activate
conda config --set auto_activate_base true
find *log |wc -l
#备份数据 【原始数据不动original】
cp -r 00data 00data2
#数据上次确认
ls -l | grep ".gz$" > 1.txt
查看虚拟环境列表
conda env list
创建虚拟环境,防污染环境变量,如果有的软件在Solving environment步骤数小时无法安装,可以新建环境
conda create -n meta
加载环境
conda activate meta
### 质量评估fastqc
# =为指定版本,-c指定安装源,均可加速安装
# -y为同意安装
conda install fastqc=0.11.9 -c bioconda -y
fastqc -v
### 评估报告汇总multiqc
# 注1.7为Python2环境,1.8/9新版本需要Python3的环境
conda install multiqc=1.9 -c bioconda -y
multiqc --version
### 质量控制流程kneaddata
conda install kneaddata=0.7.4 -c bioconda -y
kneaddata --version
trimmomatic -version # 0.39
bowtie2 --version # 2.4.2
db= /home/dengqr/dataset/metagenome/ #这里直接cd 到要放的目录吧
# 查看可用数据库
kneaddata_database
# 包括人基因组bowtie2/bmtagger、人类转录组、核糖体RNA和小鼠基因组
# 下载人基因组bowtie2索引 3.44 GB
mkdir -p $/home/dengqr/dataset/metagenome/kneaddata/human_genome
kneaddata_database --download human_genome bowtie2 $/home/dengqr/dataset/metagenome/kneaddata/human_genome
# 数据库下载慢或失败,附录有百度云和国内备份链接
下载到了 $/home/dengqr/home/dengqr/dataset/metagenome/kneaddata/
#mv /home/dengqr/$/home/dengqr/dataset/metagenome/kneaddata/ /home/dengqr/dataset/metagenome/ 【已完成】
## 1.2 (可选)FastQC质量评估
# 第一次使用软件要记录软件版本,文章方法中必须写清楚
fastqc --version # 0.11.8
# time统计运行时间,fastqc质量评估
# *.gz为原始数据,-t指定多线程
time fastqc seq/*.gz -t 2
▼ #32线程time= 27分钟 tip:在往上加
time fastqc -o 01fastqc 00data/*.gz -t 32
multiqc将fastqc的多个报告生成单个整合报告,方法批量查看和比较
# 记录软件版本
multiqc --version # 1.5
# 整理seq目录下fastqc报告,输出multiqc_report.html至result/qc目录
multiqc -d seq/ -o result/qc
▼
multiqc -d 01fastqc/ -o 02fastqc_result/
查看右侧result/qc目录中multiqc_report.html
#去除宿主
#索引目录 "/home/dengqr/dataset/metagenome/kneaddata/human_genome/hg37dec_v0.1.1.bt2"
kneaddata -h
#去宿主后双端不匹配——序列改名 先检查下
zcat 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R2_001.fastq.gz |head -n 6
zcat 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R1_001.fastq.gz |head -n 6
cp 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R2_001.fastq.gz 00datatrain\
cp 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R1_001.fastq.gz 00datatrain\
zcat 00datatrain/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R2_001.fastq.gz |head -n 6
zcat 00datatrain/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R1_001.fastq.gz |head -n 6
(可选) 序列改名,解决NCBI SRA数据双端ID重名问题,详见[《MPB:随机宏基因组测序数据质量控制和去宿主的分析流程和常见问题》](https://mp.weixin.qq.com/s/ovL4TwalqZvwx5qWb5fsYA)。
gunzip 00datatrain/*.gz
sed -i '1~4 s/$/\\1/g' 00datatrain/*R2_001.fastq
sed -i '1~4 s/$/\\2/g' 00datatrain/*R1_001.fastq
# 再次核对样本是否标签有重复
zcat 00datatrain/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R2_001.fastq |head -n 6
zcat 00datatrain/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R1_001.fastq |head -n 6
# 结果压缩节省空间
gzip seq/*.fq
# pigz是并行版的gzip,没装可使用为gzip
pigz seq/*.fq
time kneaddata -i 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R2_001.fastq.gz \
-i 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R1_001.fastq.gz \
-o 03qc -v -t 32 --remove-intermediate-output \
--reorder --bowtie2-options "--very-sensitive --dovetail" \
-db kneaddata/human_genome
### Java不匹配——重装Java运行环境
若出现错误 Unrecognized option: -d64,则安装java解决:
conda install -c cyclus java-jdk
/home/dengqr/miniconda2/bin/trimmomatic
type trimmomatic
type fastqc
"/home/dengqr/miniconda3/share/trimmomatic-0.39-2/"
time kneaddata -i 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R2_001.fastq.gz \
-i 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R1_001.fastq.gz \
-o 03qc -v -t 32 --remove-intermediate-output \
-db kneaddata/human_genome
time kneaddata -i 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R2_001.fastq.gz \
-i 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R1_001.fastq.gz \
-o 03qc -v -t 32 --remove-intermediate-output \
--trimmomatic home/dengqr/miniconda2/bin/trimmomatic \
--reorder --bowtie2-options "--very-sensitive --dovetail" \
-db kneaddata/human_genome
time kneaddata -t 40 -v \
-i 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R2_001.fastq.gz \
-i 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R1_001.fastq.gz \
-o 03qc/ \
--trimmomatic /home/dengqr/miniconda3/share/trimmomatic-0.39-2/ \
--max-memory 80g \
--trimmomatic-options "SLIDINGWINDOW:4:20 MINLEN:50" \
-db kneaddata/human_genome/ \
--bowtie2-options "--very-sensitive --dovetail --reoeder" \
--remove-intermediate-output
"/home/dengqr/miniconda2/bin/trimmomatic-0.33.jar"
"/home/dengqr/miniconda3/bin/trimmomatic"
▼
time kneaddata \
-i 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R2_001.fastq.gz \
-i 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R1_001.fastq.gz \
-o temp/qc -v -t 40 --remove-intermediate-output \
--trimmomatic /home/dengqr/miniconda3/share/trimmomatic-0.39-2/ \
--trimmomatic-options "SLIDINGWINDOW:4:20 MINLEN:50" \
--reorder --bowtie2-options "--very-sensitive --dovetail" \
-db kneaddata/human_genome
▼
time kneaddata \
-i 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R2_001.fastq.gz \
-i 00data/OSCC35A_20211015NA_AGGCAGAA_S156_L002_R1_001.fastq.gz \
-o temp1/qc -v -t 40 --remove-intermediate-output \
--trimmomatic /home/dengqr/miniconda3/share/trimmomatic-0.39-2/ \
--trimmomatic-options "SLIDINGWINDOW:4:20 MINLEN:50" \
--reorder --bowtie2-options "--very-sensitive --dovetail" \
-db kneaddata/human_genome
▼
# 采用kneaddata附属工具kneaddata_read_count_table
kneaddata_read_count_table --input temp1/qc \
--output temp1/kneaddata.txt
# 筛选重点结果列
cut -f 1,2,4,12,13 temp1/kneaddata.txt | sed 's/_1_kneaddata//' > temp1/qc/sum.txt
cat temp1/qc/sum.txt
#质控结果
fastqc temp1/qc/*_1_kneaddata_paired_*.fastq -t 2 -o temp1
multiqc -d temp1/ -o temp1/
fastqc temp1/qc/*R2_001_kneaddata_paired_*.fastq -t 2 -o temp1
multiqc -d temp1/ -o temp1/
OSCC35A_20211015NA_AGGCAGAA_S156_L002_
"/home/dengqr/dataset/metagenome/00data/Control105A_R1_001.fastq.gz"
#多任务并行运行
→记得知情同意 # 打will cite承诺引用并行软件parallel
parallel --citation
parallel -j 3 --xapply "echo 00data/{1}_R1_001.fastq.gz 00data/{1}_R2_001.fastq.gz" ::: `tail -n+2 metadata.txt|cut -f1`
time parallel -j 2 --xapply \
"kneaddata -i 00data/{1}_R1_001.fastq.gz \
-i 00data/{1}_R2_001.fastq.gz \
-o temp/qc -v -t 40 --remove-intermediate-output \
--trimmomatic /home/dengqr/miniconda3/share/trimmomatic-0.39-2/ \
--trimmomatic-options 'SLIDINGWINDOW:4:20 MINLEN:50' \
--reorder --bowtie2-options '--very-sensitive --dovetail' \
-db kneaddata/human_genome" ::: `tail -n+2 metadata.txt|cut -f1`
##
# 安装conda-pack
conda install -c conda-forge conda-pack
#conda 安装不成功可用pip安装
# pip conda-pack
# 安装好的环境下打包导出
conda pack -n humann2 -o humann2.tar.gz
# 下载
wget -c http://210.75.224.110/db/humann2/humann2.tar.gz
# 新建文件夹存放humann2环境 "/home/dengqr/miniconda3/envs/"
mkdir -p /home/dengqr/miniconda3/envs/humann2
tar -xzf humann2.tar.gz -C /home/dengqr/miniconda3/envs/humann2
# 激活环境
source /home/dengqr/miniconda3/envs/humann2/bin/activate
测试流程是否可用 time=20分钟
humann2_test
humann2 --version
切换到cd mi 3/en /hu
humann2_databases --download utility_mapping full /home/dengqr/miniconda3/envs/humann2
# 微生物泛基因组 5.37 GB 【要完整路径】【手动创文件夹】
humann2_databases --download chocophlan full /home/dengqr/miniconda3/envs/humann2/chocophlan
# 功能基因diamond索引 10.3 GB 【注意文件夹位置】
humann2_databases --download uniref uniref90_diamond /home/dengqr/miniconda3/envs/humann2
#备用下载【注意切换】 →选择百度云
https://github.com/YongxinLiu/MicrobiomeStatPlot/blob/master/Data/BigDataDownlaodList.md
拉进去/home/dengqr/miniconda3/envs/humann2
tar xvzf uniref90_annotated_1_1.tar.gz
#查看数据框是否配置上了
humann2_config
mkdir temp/concat
# 双端合并为单个文件
for i in `tail -n+2 metadata.txt|cut -f1`;do
cat temp/qc/${i}_1_kneaddata_paired_?.fastq \
> temp/concat/${i}.fq;
done
# 查看样品数量和大小
ls -sh temp/concat/*.fq