跟着Bioinformatics学数据分析:StainedGlass可视化展示基因组水平上的tandem repeat

用户7010445

发布于 2023-08-23 10:47:44

4570

发布于 2023-08-23 10:47:44

文章被收录于专栏：小明的数据分析笔记本小明的数据分析笔记本

论文

StainedGlass: interactive visualization of massive tandem repeat structures with identity heatmaps

代码链接

https://mrvollger.github.io/StainedGlass/

https://github.com/mrvollger/StainedGlass

这个工具是用来可视化展示基因组水平上tandem repeat 的相似性，是用snakemake搭建的一个流程，今天的推文我们试着拆解一下这个流程里都有哪些步骤

这个流程依赖的软件是通过搭配conda配置文件的方式去安装，但是在集群上的计算节点很多时候是不能联网的，所以最好还是提前配置好依赖软件，依赖的软件在 workflow/env目录下的env.yaml和R.yaml下

  - pandas
  - numpy
  - numba
  - cooler
  - minimap2==2.18
  - bedtools
  - samtools>=1.9
  - pysam
  - snakemake>=7.8
  - snakefmt
  - bwa
  - pigz

  - xorg-libx11
  - xorg-libxau
  - r-base>=4.0
  - r-essentials
  - r-cairo
  - r::r-tidyverse
  - r-data.table
  - r-cowplot
  - r-argparse>=2.1.2
  - r-glue
  - r::r-rcolorbrewer
  - r::r-scales
  - r::r-ggplot2
  - r-r.utils

把依赖的软件和R包都安装一下

运行命令

snakemake -s ~/biotools/StainedGlass/workflow/Snakefile --configfile=/home/myan/biotools/StainedGlass/config/config.yaml --config sample=A fasta=/data/myan/raw_data/practice/stainedGlass/chr8_cen.fasta --cores 8 make_figures -pn

会展示出这个流程每一步具体执行的命令，然后我们分别执行其中的命令看看每一步具体做了什么事

首先是对输入数据进行索引

samtools faidx chr1.fa

bedtools利用fai文件生成bed文件

## -s 参数可以设置滑窗 -w设置的是步长
bedtools makewindows -g chr1.fa.fai -w 2000 > output.bed

bedtools根据bed文件分隔fasta文件

bedtools getfasta -fi chr1.fa -bed output.bed > output.2000.fasta

batch_bed_files.py 这个脚本的作用好像是把bed文件进行分隔，--outputs参数后好像可以自定义写多少个输出

python ../batch_bed_files.py output.bed --outputs a0.bed a1.bed a2.bed

这一步是对参考构建数据库

minimap2 -f 1000 -s 400 -ax ava-ont -d output.fasta.mmi output.2000.fasta

这里的-f和-s参数没看懂是什么意思

minimap2的帮助文档

image.png

根据分隔的bed文件分别提取fasta序列

bedtools getfasta -fi chr1.fa -bed a0.bed > a0.fa
bedtools getfasta -fi chr1.fa -bed a1.bed > a1.fa
bedtools getfasta -fi chr1.fa -bed a2.bed > a2.fa

minimap2比对生成bam文件并合并

minimap2 -t 4 -f 10000 -s 400 -ax ava-ont --dual=yes --eqx output.fasta.mmi a0.fa | samtools sort -m 4G -o a0.bam
minimap2 -t 4 -f 10000 -s 400 -ax ava-ont --dual=yes --eqx output.fasta.mmi a1.fa | samtools sort -m 4G -o a1.bam
minimap2 -t 4 -f 10000 -s 400 -ax ava-ont --dual=yes --eqx output.fasta.mmi a2.fa | samtools sort -m 4G -o a2.bam

samtools merge -@ 4 -O BAM merged.bam a0.bam a1.bam a2.bam
samtools index merged.bam

接下来是画图，这里的两个python脚本起到了什么作用暂时还不太明白

python samIdentity.py --threads 8 --matches 400 --header merged.bam > output.tbl

bgzip -c output.tbl > output.tbl.gz

python refmt.py --window 2000 --fai chr1.fa.fai --full output.full.tbl.gz output.tbl.gz full.bed.gz

mkdir -p results/abc_figures/pdfs
mkdir -p results/abc_figures/pngs

Rscript aln_plot.R -b full.bed.gz --threads 8 --prefix abc

输出的部分结果