生信小白记录3-搭建上游平台踩坑记录（SNV注释）

搭建上游平台踩坑记录
test
一、质量控制
新服务器上搭建失败，conda版本过旧很多依赖包不匹配，转换宏基因组平台搭建（4.5.x，现在基本都22.x）
宏基因组平台3.12python出问题重装3.8版本，安装kneaddata，下载参考基因组hg38并解压，安装bowtie2对参考基因组构建索引bowtie2-build hg38.fa hg38
直接运行脚本中命令出问题，首先是line（循环运行需求）对不上，后续在考虑，目前先只测试一个样本，发现安装kneaddata附带的trimmomatic出问题，发现是运行命令时里面会自动指定用java打开/home/professor/anaconda3/envs/humy/bin下的trimmomatic，然而事实上该文件只是一个跳转文件，于是在脚本命令中增加了--trimmomatic /home/professor/anaconda3/envs/humy/share/trimmomatic-0.39-2指定trimmomatic所在的目录，该目录下有trimmomatic和trimmomatic.jar两个文件，这里我将trimmomatic改名为trimmomatic__，以保证运行脚本命令时启动trimmomatic时是用java命令打开trimmomatic.jar，最终运行成功。
time kneaddata --input1 rawdata/CSM7KOPI_R1.fastq.gz --input2 rawdata/CSM7KOPI_R2.fastq.gz -t 28 -db Database/hg38/Homo_sapiens/hg38 --output 01Knead/ --output-prefix CSM7KOPI --trimmomatic /home/professor/anaconda3/envs/humy/share/trimmomatic-0.39-2 --trimmomatic-options 'SLIDINGWINDOW:4:20 MINLEN:50 LEADING:3 TRAILING:3'
二、物种注释
安装kraken2（需要conda版本更新，只能在linux环境中安装，以下是在自己电脑尝试）, 构建注释数据库kraken2-build --standard --db /mnt/f/Linux/kraken2_Database
过程中出现报错Masking low-complexity regions of downloaded library...k2mask: /mnt/f/Linux/anaconda/share/kraken2-2.1.3-0/libexec/../../../lib/libstdc++.so.6: version `GLIBCXX_3.4.29' not found (required by k2mask)
通过 strings /mnt/f/Linux/anaconda/lib/libstdc++.so.6 | grep GLIBCXX命令可以检查libstdc++文件支持GLIBCXX的版本，尝试后发现
/mnt/f/Linux/anaconda/lib/libstdc++.so.6不支持3.4.29及以上，而strings /usr/lib/x86_64-linux-gnu/libstdc++.so.6 | grep GLIBCXX支持3.4.29，3.4.30
用cp命令将/usr/lib/x86_64-linux-gnu/libstdc++.so.6文件复制粘贴到了/mnt/f/Linux/anaconda/lib文件夹下（将原来的文件重命名做备份）
重新运行命令kraken2-build --standard --db /mnt/f/Linux/kraken2_Database
终止两次 原因分别是网络连接出错和占用空间满了（0.3T），于是扩展硬盘重新运行命令，它会重新下载library文件夹中的文件
library中bacteria特别大160多个G 下完这个下别的很慢 但运行上个命令会重新下载，于是换别的命令，下载archaea viral fungi human UniVec_Core对应的library
kraken2-build --download-library archaea --db /mnt/f/Linux/kraken2_Database
依次下载完这几个库后，运行以下命令建立索引
kraken2-build --build --db /mnt/f/Linux/kraken2_Database
构建中断，原因是需要很大的运行内存，Failed attempt to allocate 79063256500bytes;于是先把下载的数据库先传到服务器上再重新运行
##########################################################################################################################
发现师兄在新服务器上都下载好了 数据库 在/public/bio/database/文件夹下。
重新开始，新建环境humy
一、质量控制
安装 kneaddata，conda install -c bioconda kneaddata。
mkdir -p 01Knead
time kneaddata --input1 rawdata/CSM7KOPI_R1.fastq.gz --input2 rawdata/CSM7KOPI_R2.fastq.gz -t 28 -db /home/stu_6/HMY/hg38/ --output 01Knead/ --output-prefix CSM7KOPI --trimmomatic /home/stu_6/anaconda3/envs/humy/share/trimmomatic-0.39-2/ --trimmomatic-options 'SLIDINGWINDOW:4:20 MINLEN:50 LEADING:3 TRAILING:3'
#污染去除（先跳过了，需要相应视为污染的数据库）
kneaddata --input1 01Knead/CSM7KOPI_paired_1.fastq --input2 01Knead/CSM7KOPI_paired_2.fastq -t 28 -db Database/contaminate_db --output 01Decom/ --output-prefix CSM7KOPI --bypass-trim

二、物种注释（一、→）
# 第一步，使用Kraken2分析，不使用--use-mpa-style模式，运行失败报错invalid option -- 'H' Usage: classify [options] <fasta/fastq file(s)>。
原因是conda版本低，下载的kraken2版本是2.0.7beta_xxxxx_0，我去网站下了2.0.7beta_xxxx_3安装包自行安装(更高的版本不适配环境)，然后运行成功。
mkdir -p Kraken2_out
time kraken2 --db /public/bio/database/kraken2_230605/ --threads 56 --report Kraken2_out/CSM7KOPI.report --output Kraken2_out/CSM7KOPI.output  --paired 01Knead/CSM7KOPI_paired_1.fastq 01Knead/CSM7KOPI_paired_2.fastq
# 第二步，使用Braken校正，报错bash: bracken: command not found...，Failed to search for file: /mnt/media/BaseOS was not found。
简单来说就是bracken命令没识别，我去网站下了bracken-2.9-py38h2494328_0.tar.bz2安装包自行安装，成功运行。
mkdir -p bracken_out
time bracken -d /public/bio/database/kraken2_230605/ -i Kraken2_out/CSM7KOPI.report -o bracken_out/CSM7KOPI.S.bracken -w bracken_out/CSM7KOPI.S.bracken.report -r 150 -l S
# 第三步，将Braken的report格式转换成--use-mpa-style格式，成功。
# 用Bracken的子程序，不知道conda安装的Bracken有没有该文件，我去https://github.com/jenniferlu717/KrakenTools下了一个kreport2mpa.py。
python kreport2mpa.py -r bracken_out/CSM7KOPI.S.bracken.report -o bracken_out/CSM7KOPI.final.report

三、assembly（一、→）
安装 megahit，conda install -c bioconda megahit，使用 megahit 对原始的 paired-end 测序数据进行组装，生成 contigs。参数 -1 和 -2 指定了 paired-end 数据的路径，-m 0.9 指定了最小比对长度的阈值，-t 28 指定了线程数，--presets meta-sensitive 使用了 meta-sensitive 预设参数。
mkdir -p assem
megahit -1 01Knead/CSM7KOPI_paired_1.fastq -2 01Knead/CSM7KOPI_paired_2.fastq -m 0.9 -t 28 --presets meta-sensitive --out-prefix CSM7KOPI -o assem/CSM7KOPI_assem

四、基因预测（三、→）
安装megahit， conda install -c bioconda prodigal，使用 prodigal 对组装好的 contigs 进行基因预测，生成蛋白质和核酸序列，同时生成 GFF 文件。参数 -a 指定蛋白质序列的输出路径，-i 指定输入 contigs 文件路径，-d 指定核酸序列的输出路径，-o 指定 GFF 文件的输出路径，-f gff 指定输出格式为 GFF，-p meta 使用 meta 参数。
mkdir -p Prodigal 
prodigal -a Prodigal/CSM7KOPI_prot.faa -i assem/CSM7KOPI_assem/CSM7KOPI.contigs.fa -d Prodigal/CSM7KOPI_nucl.fa -o Prodigal/CSM7KOPI_gene.gff -f gff -p meta -s Prodigal/CSM7KOPI_stat.txt

五、生成非冗余基因集（四、→）（这步输入应该是拼接所有样本的_nucl.fa文件，这里只是测试，故输入CSM7KOPI_nucl.fa）
安装cd-hit，conda install -c bioconda cd-hit使用 cd-hit-est 对所有基因进行聚类，生成非冗余基因集。参数 -i 指定输入基因文件，-o 指定输出非冗余基因集文件，-aS 0.9 指定比对相似性的阈值，-c 0.95 指定聚类的阈值，-T 28 指定线程数。
mkdir -p cd-hit
cd-hit-est -i Prodigal/CSM7KOPI_nucl.fa -o cd-hit/genes_all_unique.fa -aS 0.9 -c 0.95 -G 0 -g 0 -T 28 -M 0

六、基因功能注释（四、→）
安装eggnog-mapper，conda install -c bioconda eggnog-mapper，找到emapper.py所在位置。使用 emapper.py 对基因进行功能注释。首先，运行一条命令对输入文件进行注释，然后使用注释结果运行另一条命令。参数 -m diamond 指定使用 Diamond 进行注释，--no_annot 不生成注释文件，--no_file_comments 不在输出文件中添加注释，-cpu 28 指定线程数，-data_dir 指定数据库目录。第一条命令报错，原因是diamond版本过低只有0.9，我去网站下了diamond-2.0.5-h43eeafb_0.tar.bz2自行安装（过高不行，不适配），
mkdir -p egg_func
/home/stu_6/anaconda3/envs/humy/bin/emapper.py -m diamond --cpu 28 --data_dir /public/bio/database/emapper_2110/ -i Prodigal/CSM7KOPI_prot.faa -o egg_func/CSM7KOPI
(/home/stu_6/anaconda3/envs/humy/bin/emapper.py --annotate_hits_table egg_func/CSM7KOPI.emapper.seed_orthologs --cpu 28 --data_dir /public/bio/database/emapper_2110/ -o egg_func/CSM7KOPI_)(这步感觉不需要 因为原本第一条命令有--no_annot，我把它去掉了，就会生成注释文件，和第二条命令产生的文件一样，所以第二条不需要运行)

七、基因丰度（三、→）
安装coverm，conda install -c bioconda coverm，使用 coverm 计算 contig 的相对丰度。参数 -t 28 指定线程数，-r 指定参考 contig 文件，--coupled 指定输入 paired-end 测序数据，--bam-file-cache-directory 指定 BAM 文件缓存目录。使用CoverM通过将reads映射到非冗余参考来估计基因丰度，并计算原始contigs中基因的覆盖率
coverm contig -t 28 -r assem/CSM7KOPI_assem/CSM7KOPI.contigs.fa --coupled 01Knead/CSM7KOPI_paired_1.fastq 01Knead/CSM7KOPI_paired_2.fastq --bam-file-cache-directory coverm_bamout
samtools view coverm_bamout/CSM7KOPI.contigs.fa.CSM7KOPI_paired_1.fastq.bam | head（查看序列比对结果）

八、snv calling（一、→，使用可以参考https://midas2.readthedocs.io/en/latest/snv_module.html）
https://github.com/czbiohub-sf/MIDAS2下载midas2包，提示需要python3.9，于是转而安装MIDAS2-1.0.2，conda install cython
，去网站下载pysam-0.21.0-py37hee149a5_0.tar.bz2自行安装，去网站https://files.pythonhosted.org/packages/b2/0a/eddd2f6667198100adb8fa6f983db66022e3c375c40d19795a3de05df986/gffutils-0.12-py3-none-any.whl#sha256=35aef6d02616ee658e846052a0d2eee8dd10a17aefc9829305fc5a28878d7702下载gffutils的whl，pip install后再次运行python setup.py install成功装包。运行midas2报错ImportError: libdeflate.so.0: cannot open shared object file: No such file or directory。我在/home/stu_6/anaconda3/envs/humy/lib/中找到libdeflate.so文件复制粘贴了一份libdeflate.so.0.
#物种选择
单样本分析
midas2 run_species --sample_name CSM7KOPI -1 01Knead/CSM7KOPI_paired_1.fastq -2 01Knead/CSM7KOPI_paired_2.fastq --midasdb_name uhgg --midasdb_dir /public/bio/database/midas2_uhgg_v110/ --num_cores 4 midas2_output
（跨样本合并，测试没有使用。midas2 merge_species --samples_list list_of_samples.tsv --min_cov 2 midas2_output/merge）
#SNV分析
单样本分析
midas2 run_snps --sample_name CSM7KOPI -1 01Knead/CSM7KOPI_paired_1.fastq -2 01Knead/CSM7KOPI_paired_2.fastq --midasdb_name uhgg --midasdb_dir /public/bio/database/midas2_uhgg_v110 --select_by median_marker_coverage,unique_fraction_covered --select_threshold=2,0.5 --num_cores 8 midas2_output
交叉样本分析
（跨样本合并计算总体SNV，测试没有使用。midas2 merge_snps --samples_list list_of_samples.tsv --midasdb_name uhgg --midasdb_dir my_midasdb_uhgg --num_cores 8 midas2_output/merge）
#拷贝数（CNV）变异分析
单样本分析
midas2 run_genes --sample_name CSM7KOPI -1 01Knead/CSM7KOPI_paired_1.fastq -2 01Knead/CSM7KOPI_paired_2.fastq --midasdb_name uhgg --midasdb_dir /public/bio/database/midas2_uhgg_v110 --select_by median_marker_coverage,unique_fraction_covered --select_threshold=2,0.5 --num_cores 8 midas2_output
（跨样本合并，测试没有使用。midas2 merge_genes --samples_list list_of_samples.tsv --midasdb_name uhgg --midasdb_dir my_midasdb_uhgg --num_cores 8 midas2_output/merge）

1360MB 质控+SNV 15:58开始 species 16:12 SNV 16:30 约35分钟 CNV 17.48 约2个小时