将测序后的reads与组装好的基因组做alignment（校准），这个过程就被叫做mapping。Mapping之后生成的SAM/BAM文件，可以获取reads mapping回参考基因组的信息（比如mapping rate，coverage，depth），从而评估基因组组装的质量。
1.Mapping工具
| reads | mapping tools |
|---|---|
| Illumina DNA-seq reads | BWA |
| Pacbio reads/ONT reads | minimap2 |
| Illumina RNA-seq | HiSat2 |
2.评估指标
（1）mapping rate
（2）Genome coverge
（3）Depth
3.开始mapping！
（1）用BWA-MEM+samtools对Illumina reads二代进行mapping
首先下载bwa（试了官网上的git clone,有点不兼容，还是选择conda）
conda install -c bioconda bwa
建立索引，使用的ref.fa是前面步骤中生成的组装好改名好的染色体fa文件，即Ctun.chr.fa。最终生成了5个文件。
bwa index ref.fa
这里的r1.clean.fq为survey数据。
bwa mem -t 32 ref.fa r1.clean.fq r2.clean.fq | samtools sort -@ 8 -m 4G > illumina.bam
或者还可以使用minimap2进行分析。两种方式二选一即可
minimap2 -t 32 -ax sr ref.fa read1.fq read2.fq | samtools sort -@ 8 -m 4G >illumina.bam
（2）用minimap2对三代reads长读长（PacBio/Nanopore reads）进行mapping
只有一种hifi数据就只用一种。
minimap2 -t 32 -ax map-hifi ../M3A_renew.genome.fa ../../rawdata/hifi1/M3-1.hifi.fastq.gz ../../rawdata/hifi2/m64292e_220707_145700.fastq.gz > M3A_hifimap.sam
我暂时没有ont数据，这行不运行
minimap2 -t 32 -ax map-ont ref.fa ont_reads.fq.gz > nanopore.sam
再生成一个bam文件
minimap2 -t 32 -ax map-hifi ../M3B_renew.genome.fa ../../rawdata/hifi1/hifi.fastq.gz ../../rawdata/hifi2/hifi.fastq.gz | samtools sort -@8 -m 4G > ..hifi.bam
（3）用HiSat2对RNA-seq数据进行mapping
create -n hisat2 hisat2
建立索引，这一步会产生8个索引文件
hisat2-build ref.fa ref.hisat
将有的数据（例如花、叶的转录组数据合并到一个bam文件中)
hisat2 --dta -p 32 -x ref.index -1 rna1_1.fa -2 rna1_2.fa 2>rna1_hisat.log | samtools
sort -@ 32 > rna1_hisat.bam
hisat2 --dta -p 32 -x ref.index -1 rna2_1.fa -2 rna2_2.fa 2>rna2_hisat.log | samtools
sort -@ 32 > rna2_hisat.bam
合并多个bam文件到一个bam文件
samtools merge -@ 32 merged_hisat.bam rna1_hisat.bam rna2_hisat.bam
4.以部分量化指标评估组装质量
（1）利用bamdst统计
bamdst安装
git clone https://github.com/shiquan/bamdst.git
cd bamdst
make
查看相关帮助
./bamdst -h
bamdst需要bed注释文件，可从fai索引文件转换
samtools faidx M3A_renew.genome.fa
cat M3A_renew.genome.fa.fai | awk '{print $1"\t1\t"$2}' > M3A.bed
run bamdst
~/pack/bamdst/bamdst -p M3A.bed M3A_hifi.bam -o bamdst.M3A_hifi.out
然而，如果基区域太大，bamdst将无法运行，内存容量无法支持。
（3）分别进行统计-1（我的办法）
A.mapping rate
samtools flagstat illumina.bam > illumina.flagstat
结果分析
674611299 + 0 in total (QC-passed reads + QC-failed reads)#共有674611299条reads通过QC+0条reads未通过QC，后面的信息行中+后的都是代表QC没通过的reads的数量。674611299 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
8019221 + 0 supplementary
0 + 0 duplicates
654899152 + 0 mapped (97.08% : N/A) #mapping record rate
666592078 + 0 paired in sequencing
333296039 + 0 read1 # 双端reads中read1的总数
333296039 + 0 read2 # 双端reads中read2的总数
573699762 + 0 properly paired (86.06% : N/A) # 86.06%比例的reads成对的映射上
644316920 + 0 with itself and mate mapped# read映射上但配对read没映射上的数量
2563011 + 0 singletons (0.38% : N/A)# 0.38%比例的read没映射上的同时，配对read映射上
62718324 + 0 with mate mapped to a different chr# reads和配对reads映射到不同染色体的情况下的reads数量
41749342 + 0 with mate mapped to a different chr (mapQ>=5)# reads和配对reads映射到不同染色体，且映射质量大于等于5的情况下的reads数量
需要注意的是，同一reads可能多次mapping，有多条记录，所以recorder number的数量会比reads number多，这里应该进行一步计算：
Mapping rate=（654899152-0）/（674611299-0）=97.08% secondary那行的数据
B.breadth of coverage
samtools depth -aa sample.bam >depth.out #计算所有单个位点的深度
cat illnmina_depth.out|awk '{x[$3]++;} END{for(i in x) print(i ":" x[i])}'#统计出深度为0的碱基位点总数量
breadth of coverage=1-深度为0碱基数/所有碱基数
C.average depth
下载mosdepth。利用mosdepth计算出的结果与老师给出的方法不一致，推测是算法不一样。
conda -n mosdepth mosdepth
生成bam文件的索引文件sample.bam.bai
samtools index sample.bam
计算深度
mosdepth -t 4 out sample.bam
输出结果为四个文件
out.mosdepth.summary.txt：文件包含每条染色体和整个基因组的信息，长度，mapped 碱基数量，平均深度，最小深度和最大深度
out.mosdepth.global.dist.txt：文件包含累积分布，指示给定覆盖率阈值下覆盖的总碱基的比例
out.per-base.bed.gz：每个碱基的输出数据
out.per-base.bed.gz.csi
（4）分别进行统计-2（老师的办法）
不知道为什么有时候-@（线程数）会出错，我的办法是直接删除。
# mean depth
samtools depth -q 4 -a file.bam | awk '{c++;s+=$3}END{print s/c}' > file.bam.depth
# breadth of coverage
samtools depth -q 4 -a file.bam | awk '{c++; if($3>0) total+=1}END{print (total/c)*100}' > file.bam.coverage
# mapping rate（我运行着出来的结果是空白）
samtools flagstat -@ 4 file.bam | awk -F "[(|%]" 'NR == 3 {print $2}' > file.bam.maprate
（5）名词解释
覆盖度（genome covarage）或者覆盖范围（breadth of coverage），参考基因组被一定深度覆盖的碱基的百分比。例如，一个基因组的90%区域在1X深度覆盖，另外仍然有70%的区域被覆盖了5X深度。
测序深度（sequencing depth）或者覆盖深度（depth of coverage），每一碱基的覆盖率是基因组碱基被测序的平均次数。 基因组的覆盖深度是通过与基因组匹配的所有短读的碱基数目除以该基因组的长度来计算的。 它通常表示为1X、2X、3X、…(1、2或3倍覆盖)。
