画测序深度分布图

1. 将测序数据比对到参考基因组。

## make ref index ##
$ minimap2 -t 10 -d index ref.fasta

## map to ref ##
$ minimap2 -t 10 -ax map-ont ref.fasta ont-sequencing-reads.fastq.gz > alin.sam

## sam2bam ##
$ samtools view -bS alin.sam > alin.bam
rm ./alin.sam

## sort ##
$ samtools sort -o ./alin.sort.bam -@ 10 alin.bam
rm ./alin.bam

## make bam index ##
$ samtools index ./alin.sort.bam -@ 6 ./alin.sort.bam.bai

2.求出每一个碱基上的深度

## depth ##
mkdir ./depth
$ samtools depth ./alin.sort.bam > ./depth/alin.sort.depth.txt

结果如下所示（其中第一列为染色体，第二列为position，第三列即为这个position的深度）

$ head alin.sort.depth.txt 
ctg000850       1       21
ctg000850       2       21
ctg000850       3       21
ctg000850       4       22
ctg000850       5       22

3.对测序深度的depth做统计

cut -f3 alin.sort.depth.txt |sort|uniq -c > alin.sort.depth.txt.count

其中第二列为depth，第一列为这个depth的频数，如第一行，即表示有231个点的深度为10.

$ head  alin.sort.depth.txt.count
   231 10
    412 11
   1032 12
   3191 13
   8499 14

4.R画图

library(ggplot2)
mytest<-read.table("alin.sort.depth.txt.count")
names(mytest)<-c("depth","count")
ggplot(mytest,aes(x=depth,y=count))+geom_point(col="red")