論文

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配置文件的方式去安裝，但是在集群上的計算節(jié)點很多時候是不能聯(lián)網(wǎng)的，所以最好還是提前配置好依賴軟件，依賴的軟件在 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

會展示出這個流程每一步具體執(zhí)行的命令，然后我們分別執(zhí)行其中的命令看看每一步具體做了什么事

首先是對輸入數(shù)據(jù)進行索引

samtools faidx chr1.fa

bedtools利用fai文件生成bed文件

## -s 參數(shù)可以設(shè)置滑窗 -w設(shè)置的是步長
bedtools makewindows -g chr1.fa.fai -w 2000 > output.bed

bedtools根據(jù)bed文件分隔fasta文件

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

batch_bed_files.py 這個腳本的作用好像是把bed文件進行分隔，--outputs參數(shù)后好像可以自定義寫多少個輸出

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

這一步是對參考構(gòu)建數(shù)據(jù)庫

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

這里的-f和-s參數(shù)沒看懂是什么意思

minimap2的幫助文檔

image.png

根據(jù)分隔的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

輸出的部分結(jié)果

image.png

image.png

這個是論文中提供的圖

image.png

推文記錄的是自己的學習筆記，很可能存在錯誤，請大家批判著看