对于通过minimap2比对的方法鉴定微生物

整体思路是： 我目前已经清楚我的测出来的菌是属于大肠杆菌，但是不太清楚具体是属于哪一种strain，于是乎我想通过将所有NCBI里面大肠杆菌的基因组序列下载然后整合成一个基因组序列，然后再用我测出来的数据和这个基因组进行比对，找到我这个菌到底与那个strain的基因组序列match程度更高

合并基因组直接用cat命令就好，这里不再赘述

直接开始使用minimap2比对

minimap2 -ax map-ont ~/project/ONT/reference/Ecoil_refeq_complete_genome/ncbi-genomes-2021-12-25/all_Escherichia_coli_complete_genome.mmi --split-prefix ./temp_sam_/ ~/project/ONT/multi_sample_basecalling_gpu_sup/0_splict_barcode/barcode01/barcode01.fastq.gz > new_barcode01_alignment.sam

如果不加--split-prefix就会出现warnning

(nanopack) qianwj@ubuntu-NF5280M5:~/project/ONT/multi_sample_basecalling_gpu_sup/5_minimap_identify/barcode01 $ minimap2 -ax map-ont ~/project/ONT/reference/Ecoil_refeq_complete_genome/ncbi-genomes-2021-12-25/all_Escherichia_coli_complete_genome.mmi ~/project/ONT/multi_sample_basecalling_gpu_sup/0_splict_barcode/barcode01/barcode01.fastq.gz > new_barcode01_alignment.sam
[WARNING] For a multi-part index, no @SQ lines will be outputted. Please use --split-prefix.
[M::main::5.410*1.00] loaded/built the index for 2342 target sequence(s)
[M::mm_mapopt_update::5.576*1.00] mid_occ = 2391
[M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 2342
[M::mm_idx_stat::5.676*1.00] distinct minimizers: 9221189 (29.13% are singletons); average occurrences: 81.657; average spacing: 5.351; total length: 4029251866

比对结束之后是一个sam文件，怎么样看对于哪种基因组比对率最高呢？
主要参考的资料：https://hc1023.github.io/2019/09/07/linux-sam-bam/#!

可以将比对之后得到的sam文件直接进行samtools查看

grep -v '^@' new_barcode01_alignment.sam | cut -f 3 | sort | uniq -c > aln_info.txt

然后再排个序

sort -t $'\t' -k 1n aln_info.txt > sorted_genome.txt

一个师兄这边给的code是

samtools view -F260 test.bam｜cut -f3｜sort ｜uniq -c ｜sort -nr -k1 ｜head

首先这个是主要做的barcode01的数据

综合这两种最后得出来的代码是：

samtools view  -F260 new_barcode01_alignment.sam | cut -f3 | sort | uniq -c | sort -rn -k1 | head

比如对于我们测的这次的barcode01的数据最后比对出来的reads的情况

(base) qianwj@ubuntu-NF5280M5:~/project/ONT/multi_sample_basecalling_gpu_sup/5_minimap_identify/barcode01 $ samtools view  -F260 new_barcode01_alignment.sam | cut -f3 | sort | uniq -c | sort -rn -k1 | head
   1376 NZ_CP080246.1
   1275 NZ_CP035486.1
   1203 NZ_CP022686.1
    828 NZ_CP007799.1
    565 NZ_CP027340.1
    458 NZ_CP037449.1
    433 NZ_CP081010.1
    429 NZ_CP032085.1
    425 NZ_CP061206.1
    408 NZ_CP070041.1

以及对于我们这次的比对的barcode03的数据最后的reads分布情况

(base) qianwj@ubuntu-NF5280M5:~/project/ONT/multi_sample_basecalling_gpu_sup/5_minimap_identify/barcode03 $ head identify_info.txt
   1734 NZ_CP027340.1
   1360 NZ_CP025520.1
   1288 NZ_CP053607.1
   1274 NZ_CP047127.1
   1247 NZ_CP045741.1
   1160 NZ_CP017100.1
    734 NZ_CP080620.1
    710 NZ_CP081007.1
    631 NZ_CP081008.1
    628 NC_017625.1

接下来以此类推barcode02以及barcode03

awk按列求和

awk '{sum += $1};END {print sum}' test

如果要下载NCBI里面的数据的话

可以用这个链接：
https://ftp.ncbi.nlm.nih.gov/genomes/