mamba activate ncbi_genome_download

## 根据登录号下载
ncbi-genome-download  --assembly-accessions refseq.txt plant  --section refseq  --formats  fasta,gff  --progress-bar  --parallel 3   -o data1  --flat-output

## 根据物种名
ncbi-genome-download --genera "Pseudomonas reinekei" bacteria                 --section genbank  --formats  fasta,gff  --progress-bar  --parallel 3 -o genbank
ncbi-genome-download --genera "Pseudomonas reinekei" bacteria                 --section refseq  --formats  fasta,gff  --progress-bar  --parallel 3 -o refseq
## 根据属名
ncbi-genome-download --genera Pseudomonas bacteria
## 物种名或属名可以放txt中
ncbi-genome-download --genera ref.txt bacteria

比对nt库 -- 污染比对

1.比对nt库
perl /share/nas6/zhouxy/functional_modules/pollution_nt_blast/pollution_nt_blast_pip_v3.pl -fqdir 路径/2.clean_data -od 输出目录

2.占比可视化
python3 /share/nas1/yuj/script/ddradANDreseqANDgenome/contamination_visualizer.py   -i huasheng8.Simplify.Pollution.stats.xls

bwa2pip

1.0 比对参考基因组（2次比对）

1.进行第一次比对

2.将第一步比对输出为reads
for i in `ls step1/mapping_dir`;do echo  $i;mkdir reads_data2/$i;done
for i in `ls step1/mapping_dir`;do echo  $i; samtools view -b -f 0x2 -F 12  step1/mapping_dir/$i/target_genome/$i.sam  > step1/mapping_dir/$i/target_genome/$i.mapped_pairs.bam && samtools fastq -1 reads_data2/$i/${i}_L007_R1_001.fq -2 reads_data2/$i/${i}_L007_R2_001.fq   step1/mapping_dir/$i/target_genome/$i.mapped_pairs.bam   &  done

3.进行第二次比对
python3 /share/nas1/yuj/pipline/find_genome/bwa2pip-v3.py  -i  reads_data2  -g  genome_data2   -o  step2   --keep_sam  --build_index 
## 如果已经建好索引，可以直接
python3 /share/nas1/yuj/pipline/find_genome/bwa2pip-v3.py  -i  reads_data2  -g  genome_data2   -o  step2   --keep_sam

4.
cat step2/map.stat.xls | sort  -k1,1 -k3,3nr

1.1 共用参考（可能不能用）

0.创建文件夹
data_path=分发路径
samples=`ls $data_path/2.clean_data/  |grep -v "md5.txt"` && echo $samples
for i in $samples;do echo $i;touch $i.txt;mkdir $i;done
mkdir summary/data/subset -p

或者
for i in $samples;do echo $i;mkdir $i;done

1.截取250万条reads（单端125万，对应500_0000行）
# zcat   xx路径/13_S52_L002_R1_001.fastq.gz |head -n 5000000 > summary/data/subset/13_S52_L002_R1_001.fq
## for循环   公司测序格式
for i in $samples;do    for j in {1,2} ; do  zcat $data_path/2.clean_data/$i/${i}_S*_L00*_R${j}_00*.fastq.gz  |head -n 5000000 > summary/data/subset/${i}_L007_R${j}_001.fq  & done; done

或者
## for循环   其他公司测序格式
for i in $samples;do    for j in {1,2} ; do  zcat $data_path/2.clean_data/$i/${i}_R${j}.fq.gz  |head -n 5000000 > $i/${i}_L007_R${j}_001.fq  & done; done

2.准备基因组文件

3.开始比对
项目目录】
for i in $samples;do echo $i;  perl /share/nas1/yuj/pipline/find_genome/bwa2pip.pl -i  summary/data/subset  -g  genome_data -o $i.ident -p $i & done

4.查看结果
for i in $samples;do head -n 1 $i.ident/map.stat.xls;done | sort -k 3 -n > ref_sort.txt
cat ref_sort.txt

1.2 不同样品对应不同参考

0.创建文件夹
data_path=分发路径
samples=`ls $data_path/2.clean_data/  |grep -v "md5.txt"` && echo $samples
for i in $samples;do echo $i;touch $i.txt;mkdir $i;done
mkdir summary/data/subset -p

1.截取250万条reads
# zcat   xx路径/13_S52_L002_R1_001.fastq.gz |head -n 5000000 > summary/data/subset/13_S52_L002_R1_001.fq
## for循环
for i in $samples;do    for j in {1,2} ; do  zcat $data_path/2.clean_data/$i/${i}_S*_L00*_R${j}_00*.fastq.gz  |head -n 5000000 > summary/data/subset/${i}_L007_R${j}_001.fq  & done; done

2.下载参考基因组
## 把样品对应的属名物种名放进样品txt中
mamba activate ncbi_genome_download
## for循环   每个文件夹下有对应txt
for i in {UB002};do mv $i.txt  $i; done

--- refseq数据库
for i in {UB002};do cd $i && ncbi-genome-download --genera $i.txt bacteria   --formats  fasta,gff  --progress-bar  --parallel 3 & done
for i in P{10P11,16P17,5}; do  cd $i && ncbi-genome-download --genera *.txt fungi   --formats  fasta,gff  --progress-bar  --parallel 3 & done
--- genbank数据库
for i in P{10P11,16P17,5}; do cd $i  && ncbi-genome-download --genera *.txt fungi   --formats  fasta,gff  --progress-bar  --parallel 3  --section genbank  -o ./ & done 

3.比对运行
cd UB002 && perl /share/nas2/yuj/project/2024/re-sequencing/GP-20240220-7911_20240509/pipline/bwa2pip.pl -i 项目路径/summary/data/subset -g refseq/bacteria -o ident -p UB002 && cd - &

## for循环
for i in {UB0029};do rm $i/ident -rf ;done
for i in {UB0029};do echo $i; cd  $i && perl /share/nas2/yuj/project/2024/re-sequencing/GP-20240220-7911_20240509/pipline/bwa2pip.pl -i 项目路径/summary/data/subset -g */fungi -o ident -p $i  && cd - & done

for i in {UB0029};do echo $i; cd  *$i* && perl /share/nas2/yuj/project/2024/re-sequencing/GP-20240220-7911_20240509/pipline/bwa2pip.pl -i 项目路径/summary/data/subset -g */fungi -o  $i.ident -p $i  && cd - & done # 不同样品的参考放一块了，需要分别对样品进行比对
# summary/data/subset里面是每个样品截取后的序列
# refseq/bacteria文件夹下有很多以登录号命名的文件夹
# */fungi文件夹下有很多以登录号命名的文件夹

4.查看结果
for i in {UB002};do head -n 1 $i/ident/map.stat.xls;done | sort -k 3 -n > ref_sort.txt
cat ref_sort.txt

## 看看有没有gff
后n个都是比对大于70%的
n=11 && for i in $(cat ref_sort.txt | tail -n $n | cut -f 2);do ll */refseq/bacteria/$i/*.gff*;done

5.没有参考基因组的，尝试genbank数据库
mamba activate ncbi_genome_download
n=11 && for i in $(cat ref_sort.txt | head -n $n | cut -f 1);do cd $i && ncbi-genome-download --genera $i.txt bacteria  --section genbank  --formats  fasta,gff --progress-bar  --parallel 3 -o ./ & done

for i in $(cat ref_sort.txt | head -n $n | cut -f 1);do echo $i; cd $i && perl /share/nas2/yuj/project/2024/re-sequencing/GP-20240220-7911_20240509/pipline/bwa2pip.pl -i 项目路径/summary/data/subset -g genbank/genbank/bacteria -o ident2 -p $i && cd - & done

for i in $(cat ref_sort.txt | head -n $n | cut -f 1);do head -n 1 $i/ident2/map.stat.xls;done | sort -k 3 -n > ref_sort2.txt
cat ref_sort2.txt

n=11 && for i in $(cat ref_sort2.txt | tail -n $n | cut -f 2);do ll */genbank/genbank/bacteria/$i/*.gff*;done # 后n个是比对大于70%的

while read -r line; do   echo "$line"; cd $(echo "$line" | cut -f 1) && ln -s refseq/bacteria/$(echo "$line" | cut -f 2) data;cd -; done < 11ok.txt

 for i in $(cat 11ok.txt | cut -f 1);do cd $i && gunzip data/*.gz && cd data && ir.py -i *.fna -l | sort -k 1 -nr > len.txt && cd - & done

 for i in $(cat 11ok.txt | cut -f 1);do cat  $i/data/len.txt ;done

方法描述

共进行了5轮基因组下载比对

1.（一轮比对）根据第1次菌种鉴定结果，
对每个样品利用ncbi-genome-download v0.3.3(https://github.com/kblin/ncbi-genome-download)从refseq数据库下载对应鉴定物种们的基因组，
从测序数据中抽取125万条reads后使用bwa-mem2 v2.2.1(https://github.com/bwa-mem2/bwa-mem2)将reads比对参考基因组，
对获得的sam文件使用samtools v1.9(https://github.com/samtools/samtools/)软件的flagstat功能进行统计获得匹配率，
从匹配率大于70％的众多基因组中挑选比对率最高的基因组作为参考。

2.（二轮比对）若某个样品没有匹配率大于70％的基因组，
再使用ncbi-genome-download从genbank数据库下载对应物种们的基因组，
重新抽取reads比对统计，从匹配率大于70％的众多基因组中挑选比对率最高的基因组作为参考。

3.若某个样品依旧没有大于70的参考，
则对该样品使用spades v3.15.5(https://github.com/ablab/spades)进行初步组装,
从组装结果中挑选100k左右的contig到ncbi使用blast进行比对获得潜在的鉴定物种，取前5个物种。

4.（三轮比对）将第3步的潜在物种们和第2次提供的菌种鉴定结果整合，
重复1步骤

5.（四轮比对）若某个样品没有匹配率大于70％的基因组，
重复2步骤

6.（五轮比对）若在第4、第5步中对某个样品操作时，遇到基因组多到软件下载不下来中途报错的情况，
则根据第4步的整合结果，将该样品对应的潜在物种们放宽至属水平，把其他样品下载下来的同属基因组整合后，
进行比对统计，挑选大于70％的参考基因组

7.若依旧没有，就没有其他操作，定为找不到