bioconvert计算测序覆盖度并使用R语言ggplot2画折

参考基因组下载自NCBI

https://www.ncbi.nlm.nih.gov/nuccore/FN433596

下载原始测序数据

最近发现了两个新方法

• 一个是bioconvert可以直接下载
• 还有一个工具是 kingfisher 这个好强大 git主页链接 https://github.com/wwood/kingfisher-download

我试了一下这个kingfisher这个工具

kingfisher get -r ERR043367 -m ena-ascp prefetch aws-http

需要提前安装Aspera这个工具并添加到环境变量

把一个软件临时添加到环境变量可以使用如下命令

先使用cd命令进入软件的可执行文件的目录

接下来是export PATH=$PWD:$PATH

计算覆盖度

首选是参考基因组构建索引

bwa index Staphylococcus_aureus.fasta

比对

bwa mem -M -t 16 Staphylococcus_aureus.fasta ERR043371_1.fastq ERR043371_2.fastq > output.sam

sam文件转换为bam

samtools view -S output.sam -O bam -o output.bam

bam文件排序

samtools sort output.bam -@ 16 -O bam -o output.sorted.bam

计算覆盖度

bioconvert bam2cov output.sorted.bam output.cov

输出文件的部分

image.png

最后是折线图可视化

df<-read.delim("output.cov",header=F)
dim(df)
head(df)
library(ggplot2)
pdf(file = "cov.pdf",width=10,height = 4)
ggplot(data=df,aes(x=V2,y=V3))+
  geom_line()+
  scale_y_continuous(expand=c(0,0))
dev.off()

image.png

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！