SparKでBAMファイルのカバレッジを見る

RNA-seq解析やDNA sequenceの結果を見るときに、よくIGV (Integrative Genomics Viewer)が使われる。

IGVのカバレッジトラックを見ることで、それぞれの領域にどのくらいのreadがmappingされたかを確認することができる。

そのカバレッジトラックのみを視覚化することのできるSparKの動かし方を以下にまとめる。

以下のサイトを参考にした。

https://kazumaxneo.hatenablog.com/entry/2020/02/05/132139

https://github.com/harbourlab/SparK


### ソフトのインストール

はじめに、condaでdeeptoolsをインストールし、さらにgithubからSparKを入手した。

```
$ conda install -c bioconda deeptools
$ git clone https://github.com/harbourlab/SparK.git
```

### データセット
<i>Aspergillus oryzae</i>のRNA-seq data (Acc. No. DRX196099, PMID: 32518660)をRIB40株のReference配列にmappingしたものをplotの描画に用いた。
以下のようにmappingし、bamファイルを用意した。
```
$ fastp -v
fastp 0.23.4
$ fastp --detect_adapter_for_pe -i DRR205700_1.fastq -I DRR205700_2.fastq -o DRR205700_clean_1.fastq -O DRR205700_clean_2.fastq -w 6

$ bowtie2 --version
version 2.5.0
$ bowtie2-build -f A_oryzae_RIB40_current_chromosomes.fasta RIB40
$ bowtie2 -x RIB40 -p 6 -1 DRR205700_clean_1.fastq -2 DRR205700_clean_2.fastq -S DRR205700.sam

$ samtools version
samtools 1.16.1
$ samtools sort -@ 6 -O bam -o DRR205700.bam DRR205700.sam 
$ samtools index DRR205700.bam 
```
### plotの描画
上記のbamファイルを用いて、カバレッジを可視化する。
bedファイルを入力することで、遺伝子領域も可視化できる。
-pr optionで出力するregionを指定する。染色体全体を可視化する場合も、数字の入力は必須。

```
$ bamCoverage -b DRR205700.bam -o output.bdg -bs 1 -of bedgraph   
            
$ python3 /Users/username/Downloads/SparK/SparK.py -pr Chr1_A_oryzae_RIB40:100000-200000 -bed transcripts.bed -cf output.bdg
```

plotは以下のようにsvgファイルとして、出力される。
![Fig.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1060100/b7fb8327-1473-3488-7d41-7e8fd57a177f.png)


#### SparKのoption

```
$ python3 /Users/username/Downloads/SparK/SparK.py -h 
usage: SparK.py [-h] [-pt PLOT_TYPE] [-ps SHOW_PLOTS] -pr REGION -cf            
                CONTROL_FILES [CONTROL_FILES ...]                               
                [-tf TREAT_FILES [TREAT_FILES ...]]                             
                [-cg CONTROL_GROUPS [CONTROL_GROUPS ...]]                       
                [-tg TREAT_GROUPS [TREAT_GROUPS ...]]                           
                [-gl GROUP_LABELS [GROUP_LABELS ...]] [-l LABELS [LABELS ...]]  
                [-gs GROUP_AUTOSCALE]
                [-es EXCLUDE_FROM_GROUP_AUTOSCALE [EXCLUDE_FROM_GROUP_AUTOSCALE ...]]
                [-eg EXCLUDE_GROUPS [EXCLUDE_GROUPS ...]]
                [-f FILLS [FILLS ...]]
                [-cs CUSTOM_Y_AXIS_SCALES [CUSTOM_Y_AXIS_SCALES ...]]
                [-dc DISPLAY_CHROM_LOCATION] [-gtf GTFFILE] [-tss DRAWTSS]
                [-genestart DRAW_GENESTART] [-sp SPARK]
                [-sc SPARK_COLOR [SPARK_COLOR ...]] [-sm SMOOTHEN]
                [-o OUTPUT_NAME] [-bed BED_FILES [BED_FILES ...]]
                [-bedcol BED_COLOR [BED_COLOR ...]]
                [-bedlab BED_LABELS [BED_LABELS ...]] [-w TRACK_WIDTH]
                [-dg DISPLAY_GENES [DISPLAY_GENES ...]]
                [-dt DISPLAY_TRANSCRIPTS [DISPLAY_TRANSCRIPTS ...]]
                [-wg WRITE_GENENAMES] [-scale DISPLAY_SCALEBAR]

SparC_args

options:
  -h, --help            show this help message and exit
  -pt PLOT_TYPE, --plot_type PLOT_TYPE
                        choices: standard, STD, sine
  -ps SHOW_PLOTS, --show_plots SHOW_PLOTS
                        choices: all, averages
  -pr REGION, --region REGION
                        example: chr1:1647389-272634
  -cf CONTROL_FILES [CONTROL_FILES ...], --control_files CONTROL_FILES [CONTROL_FILES ...]
                        separate by space
  -tf TREAT_FILES [TREAT_FILES ...], --treat_files TREAT_FILES [TREAT_FILES ...]
                        separate by space
  -cg CONTROL_GROUPS [CONTROL_GROUPS ...], --control_groups CONTROL_GROUPS [CONTROL_GROUPS ...]
                        group numbers separate by spacse
  -tg TREAT_GROUPS [TREAT_GROUPS ...], --treat_groups TREAT_GROUPS [TREAT_GROUPS ...]
                        group numbers separate by space
  -gl GROUP_LABELS [GROUP_LABELS ...], --group_labels GROUP_LABELS [GROUP_LABELS ...]
                        set group labels
  -l LABELS [LABELS ...], --labels LABELS [LABELS ...]
                        set labels for controls and treatment
  -gs GROUP_AUTOSCALE, --group_autoscale GROUP_AUTOSCALE
                        set to "yes" to autoscale all groups(plots), except
                        the groups excluded in -eg
  -es EXCLUDE_FROM_GROUP_AUTOSCALE [EXCLUDE_FROM_GROUP_AUTOSCALE ...], --exclude_from_group_autoscale EXCLUDE_FROM_GROUP_AUTOSCALE [EXCLUDE_FROM_GROUP_AUTOSCALE ...]
                        group numbers of groups to be excluded from autoscale
  -eg EXCLUDE_GROUPS [EXCLUDE_GROUPS ...], --exclude_groups EXCLUDE_GROUPS [EXCLUDE_GROUPS ...]
                        Exclude groups from the analysis
  -f FILLS [FILLS ...], --fills FILLS [FILLS ...]
                        track colors. One or two colors in hex format for
                        control and treatment tracks
  -cs CUSTOM_Y_AXIS_SCALES [CUSTOM_Y_AXIS_SCALES ...], --custom_y_axis_scales CUSTOM_Y_AXIS_SCALES [CUSTOM_Y_AXIS_SCALES ...]
                        Enter y-axis values for all groups(plots). All groups
                        need a value. Enter "D" for each group no value should
                        be assigned, e.g. to keep autoscaling functionality
                        for some groups
  -dc DISPLAY_CHROM_LOCATION, --display_chrom_location DISPLAY_CHROM_LOCATION
                        set to "no" if you do not want to plot the chromosomal
                        coordinates
  -gtf GTFFILE, --gtffile GTFFILE
                        link gtf file for drawing genes here
  -tss DRAWTSS, --drawtss DRAWTSS
                        set to "no" if TSS sites should not be indicated
  -genestart DRAW_GENESTART, --draw_genestart DRAW_GENESTART
                        set to "yes" if TSS sites should be indicated
  -sp SPARK, --spark SPARK
                        display significant change "yes"
  -sc SPARK_COLOR [SPARK_COLOR ...], --spark_color SPARK_COLOR [SPARK_COLOR ...]
                        spark color
  -sm SMOOTHEN, --smoothen SMOOTHEN
                        smoothen tracks, int
  -o OUTPUT_NAME, --output_name OUTPUT_NAME
                        output graph name, str
  -bed BED_FILES [BED_FILES ...], --bed_files BED_FILES [BED_FILES ...]
                        bed files to be plotted
  -bedcol BED_COLOR [BED_COLOR ...], --bed_color BED_COLOR [BED_COLOR ...]
                        colors of bed files in hex
  -bedlab BED_LABELS [BED_LABELS ...], --bed_labels BED_LABELS [BED_LABELS ...]
                        set labels for bed tracks
  -w TRACK_WIDTH, --track_width TRACK_WIDTH
                        width of the track, default = 150, int
  -dg DISPLAY_GENES [DISPLAY_GENES ...], --display_genes DISPLAY_GENES [DISPLAY_GENES ...]
                        genes to display from the gtf file
  -dt DISPLAY_TRANSCRIPTS [DISPLAY_TRANSCRIPTS ...], --display_transcripts DISPLAY_TRANSCRIPTS [DISPLAY_TRANSCRIPTS ...]
                        display custom transcripts. By default, all
                        transcripts annotated in the gtf file will be merged
                        and displayed as one gene. Alternatively all can be
                        plotted seperatelly by setting this to "all". Further,
                        Transcript IDs can be listed to plot only certain
                        transcripts
  -wg WRITE_GENENAMES, --write_genenames WRITE_GENENAMES
                        write genename instead of transcript ID when
                        transcripts are plotted. Set to "yes".
  -scale DISPLAY_SCALEBAR, --display_scalebar DISPLAY_SCALEBAR
                        set to "no" to remove scalebar
```