查找参考基因组

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

0.创建文件夹
mkdir summary/data/subset -p
for i in {UB002};do echo $i;touch $i.txt;mkdir $i;done

1.截取250万条reads
zcat   xx路径/13_S52_L002_R1_001.fastq.gz |head -5000000 > summary/data/subset/13_S52_L002_R1_001.fq
## for循环
for i in {UB002};do    for j in {1,2} ; do  zcat 分发路径/2.clean_data/$i/${i}_S*_L007_R${j}_001.fastq.gz  |head -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.若依旧没有，就没有其他操作，定为找不到