测序深度图的制作

1.安装samtools并将sam文件转为bam文件

samtools是一款非常好用的工具，本文我们使用的功能较少。详细的使用方法请见：samtools用法详解_samtools 使用_sunchengquan的博客-CSDN博客

conda create -n samtools
#创建安装环境
conda activate samtools
#激活环境
conda install samtools
#通过conda安装samtools
cd /to/your/work/place
#转到存放有sam文件的文件夹，sam文件为getorganelle组装叶绿体基因组时产生，位于seed文件夹下
samtools view -@8 -b test.sam > test.bam
#使用samtools view将sam文件转为bam文件,-@表示使用的线程数

2.排序归档和提取

# 构建索引
samtools sort -o test_sort.bam test.bam
#对bam文件进行排序
samtools index test_sort.bam test_sort.bam.bai
#对test_sort.bam构建索引
samtools depth test_sort.bam > test_depth.txt
#计算每一个位点或者区域的测序深度并生成文件储存
cut -f 2-3 test_depth.txt > depth.txt
#对生成的txt文件进行数据裁剪，直接使用Linux的cut命令，-f表示裁剪整列，2-3表示裁剪后保存2-3列

3.合并步长，精简信息

import os
# 脚本只需要修改file_path 和 output_file

# output_file是2000步长的合并文件
#file_path是上一步生成的深度文件depth.txt
file_path = open(r'D:/python/depth.txt', 'r')
output_file = r'D:/python/depth_2000.txt'

cont = file_path.readlines()

with open(output_file, 'w') as ff: # output file
    for i in range(len(cont)):
        if i % 2000 == 0:
            start = i
            end = min(i + 2000, len(cont))  # 确保不超过列表的长度
            middle = start + (end - start) // 2  # 计算实际区间中点
            all_sum = 0
            for j in range(start, end):
                all_sum += int(cont[j].split('\t')[2])
            average_depth = round(all_sum / (end - start), 0)
            ff.write(f"{middle}\t{average_depth}\n")

4.画测序深度图

先将上一步得到的count文件放入R的工作目录

#查看工作目录
getwd()
#设置工作目录
setwd("/path/to/your/work/place")
#如果没有安装ggplot2，可以安装install.packages(c("package1","package2",....)),此命令可同时安装多个包
install.packages("ggplot2")
#加载包
library(ggplot2)

#加载数据
data <- read.delim("D:/python/depth_2000.txt", header=FALSE, sep = '\t')
data2 <- read.delim("D:/python/depth.txt", header=FALSE, sep = '\t')

##按照2000的步长间隔进行绘图
ggplot(data, aes(V1, V2)) + 
  geom_bar(stat = 'identity', width = 800, fill = "#BBE7FC") + 
  ylim(0,1500) + theme_classic() + xlab("Sequence length") + 
  ylab("Mean base depth")+
  geom_rect(aes(xmin=1000, xmax=86319, ymin=0, ymax=0.5), 
            fill="#BBE7FC", colour="#FDDBB2", linewidth=1.5) + #xmax=86319是LSC的终点
  geom_rect(aes(xmin=86320, xmax=111890, ymin=0, ymax=1), 
            fill="#B45218", colour="#B66300", linewidth=1.5) + #xmin=86320是IRB的起点, xmax=111890是IRB的终点
  geom_rect(aes(xmin=111891, xmax=130299, ymin=0, ymax=0.5), 
            fill="#F3D68A", colour="#F4D68A", linewidth=1.5) + #xmin=111891是SSC的起点, xmax=130299是SSC的终点
  geom_rect(aes(xmin=130300, xmax=155870, ymin=0, ymax=0.5), 
            fill="#C69D3C", colour="#C69B3C", linewidth=1.5) #xmin=130300是IRA的起点, xmax=155870是IRA的终点


##按照所有位点绘图
ggplot(data2, aes(V1, V2)) + 
  geom_bar(stat = 'identity', width = 800, fill = "#BBE7FC") + 
  ylim(0,1500) + theme_classic() + xlab("Sequence length") + 
  ylab("Mean base depth")+
  geom_rect(aes(xmin=1, xmax=86319, ymin=0, ymax=0.5), 
            fill="#BBE7FC", colour="#FDDBB2", linewidth=1.5) + #xmax=86319是LSC的终点
  geom_rect(aes(xmin=86320, xmax=111890, ymin=0, ymax=1), 
            fill="#B45218", colour="#B66300", linewidth=1.5) + #xmin=86320是IRB的起点, xmax=111890是IRB的终点
  geom_rect(aes(xmin=111891, xmax=130299, ymin=0, ymax=0.5), 
            fill="#F3D68A", colour="#F4D68A", linewidth=1.5) + #xmin=111891是SSC的起点, xmax=130299是SSC的终点
  geom_rect(aes(xmin=130300, xmax=155870, ymin=0, ymax=0.5), 
            fill="#C69D3C", colour="#C69B3C", linewidth=1.5) #xmin=130300是IRA的起点, xmax=155870是IRA的终点

Rstudio可以选择输出文件的类型，有png和pdf两种。

完
本文参考链接为https://blog.csdn.net/salty_fish_xu/article/details/136592818?spm=1001.2014.3001.5502