Genome-wide off-target evaluation of adenine base editing by tracing intermediate dI.
- 2025-10-09 first commit. By Haowei Meng and Sihan Zhang.
Adenine base editors (ABEs) are CRISPR-Cas-derived tools that enable precise conversion of adenine (A) to guanine (G) in DNA without inducing double-stranded breaks. ABEs mediate A-to-G conversions through deamination of adenine to inosine (dI) as an intermediate, making the capture of dI a effective strategy to identify their off-target sites.
dI-profiling aims to find all reliable off-target edits of ABEs in the whole genome.
This set of experiments and analytical procedures were reported by the YiLab @ Peking University in 2025.
Dipro tool can help find off-target editing sites, sgRNA alignments, and performing visualization of results. One can analyze dI-profiling sequencing data with the dipro tool.
git clone https://github.com/menghaowei/dipro.git
cd dipro
# this conda env was tested on linux
conda env create -f conda-env.yaml
conda activate dipro- 1. Map the dI-profiling reads to BAM
- 2. Tracing A-to-G signal regions
- 3. Perform statistical test for candidate regions
- 4. Criteria selection and sgRNA binding site identification
All scripts can be found at /src/dipro.
This step processes raw FASTQ files into bam format, comprising following steps:
Step 1. Removing sequencing adapter in raw FASTQ
Step 2. Mapping dI-profiling reads to reference genome using BWA MEM
Step 3. Removing duplication in bam files
Step 4. Filtering reads in bam that are not properly aligned
This step catches query mutation information from bam files (i.e., A-to-G mutation signals for the Watson strand and T-to-C mutation signals for the Crick strand).
Samtools mpileup is used to generate a pileup-format output that summarizes reads alignments at each genomic position. This includes information such as base counts, quality scores and indel data.
samtools mpileup \
--reference YOUR_PATH_TO/hg38_only_chromosome.fa \
-q 20 -Q 20 \
-o mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.mpileup \
bam.bwa/dIprofiling-ABE_treat_rep1_bwa_hg38_sort_rmdup.MAPQ20.bam &An output example for the output mpileup:
chr1 3080744 G 22 ,,..,,.,,.,,.,.,.,.,., HGHGIHHHFHHFIoIHHFIEHH
chr1 3080745 G 22 ,,..,,.,,.,,.,.,.,.,., FCFFFHGHFHHHFnHFHIEFGF
chr1 3080746 A 22 ,,..,,.,,.,,.,.,.,.,., CBEEDDECC9CCEgDDFCECFD
chr1 3080747 A 22 ,,..,,.,,.,,.,.,.,.,., DCFDCDFDCDCDFgDDFDEDED
chr1 3080748 C 22 ,$,..,,.,,.,,.,.,.,.,., ECGFHGGIGHHHHpHIHGHIHH
chr1 3080749 A 21 ,..,,.,,.,,.,.,.,.,., CDFCDFCCDDCEhDDEDECED
chr1 3080750 T 21 ,..,,.,,.,,.,.,.,.,., BDCDDDDDCFDDhCECFDDDE
chr1 3080751 G 21 ,..,,.,,.,,.,.,.,.,., FHHHHIHHHHFHpHHIHHFGG
chr1 3080752 T 21 ,..,,.,,.,,.,.,.,.,., FCDDEBFDBEFCdCDCDDDBD
chr1 3080753 T 21 ,..,,.,,.,,.c.,.,.,., BCDFDCFFCFEDfDFDDDFDF
...
The parse-mpileup-ABE-v1.0.py script parses the A, G, C and T counts, as well as the indel counts, at each genomic position from mpileup, and outputs the results in bmat format.
python parse-mpileup-ABE-v1.0.py \
-i mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.mpileup \
-o mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.bmat \
-p 24 -n 0 \
> mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.bmat.log 2>&1 &Column explanation for bmat:
chr_name: str, chromosome name of signal region, value like 'chr1', 'chr2', etc.site_index: int, the coordinate of the site, and the coordinate index is based on a 1-based scale.ref_base: str, the reference base at the position, value like 'A', 'T', 'C' or 'G'.A: int, the count of reads supporting an 'A' base at the position.G: int, the count of reads supporting an 'G' base at the position.C: int, the count of reads supporting an 'C' base at the position.T: int, the count of reads supporting an 'T' base at the position.del_count: int, the number of reads with a deletion at the position.insert_count: int, the number of reads with a insersion at the position.ambiguous_count: int, the number of reads with ambiguous bases (e.g., N).deletion: str, the specific deleted sequence, e.g., 'AT' for a deletion of adenine and thymine.insertion: str, the specific inserted sequence, e.g., 'G' for an insertion of guanine.ambiguous: str, the ambiguous sequence or notation (e.g., 'N' or 'R').mut_num: int, the total number of distinct mutation types observed at the position.
An output example for the output bmat:
chr1 2529156 A 454 20 0 0 0 0 0 . . . 20
chr1 2529157 C 0 0 474 0 0 0 0 . . . 0
chr1 2529158 A 445 29 0 0 0 0 0 . . . 29
chr1 2529159 A 448 24 0 0 0 0 0 . . . 24
chr1 2529160 A 446 28 0 0 0 0 0 . . . 28
chr1 2529161 T 0 0 1 468 0 0 0 . . . 1
chr1 2529162 A 359 112 0 0 0 0 0 . . . 112
chr1 2529163 T 0 0 0 465 0 0 0 . . . 0
chr1 2529164 C 0 0 471 0 0 0 0 . . . 0
chr1 2529165 A 447 25 0 0 0 0 0 . . . 25
...
This step involves filtering query base types (i.e., A or T) from the bmat, calculating coverage, and determining the count and ratio of A-to-G or T-to-C mutation reads.
# A on the Watson strand
cat mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.bmat | awk '$3 == "A" {{print $0}}' > mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20_A.bmat
python bmat2pmat-ABE-v1.0.py \
-i mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20_A.bmat \
-o mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20_A.pmat \
--InHeader False --InLikeBED False --OutHeader True &
# T on the Watson strand
cat mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.bmat | awk '$3 == "T" {{print $0}}' > mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20_T.bmat
python bmat2pmat-ABE-v1.0.py \
-i mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20_T.bmat \
-o mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20_T.pmat \
--InHeader False --InLikeBED False --OutHeader True &Column explanation for pmat:
chr_name: str, chromosome name of signal region, value like 'chr1', 'chr2' ...site_index_start: int, the coordinate of selected site, and the coordinate index is based on a 1-based scale.site_index_end: int, the coordinate of selected site, which is exactly the same as the site_index_start value.A: int, the count of reads supporting an 'A' base at the position.G: int, the count of reads supporting an 'G' base at the position.C: int, the count of reads supporting an 'C' base at the position.T: int, the count of reads supporting an 'T' base at the position.mut_type: str, the type of dominant mutation direction at the position.ref_base: str, the reference base at the position, value like 'A', 'T', 'C' or 'G'.mut_base: str, the dominant mutation base at the position.ref_num: int, the number of reads supporting the reference base at the position.mut_num: int, the total number of reads supportingref_base.cover_num: int, the reads covering the position, calculated asmut_num+ref_num.mut_ratio: float, the ratio of mutated reads to total covering reads, calculated asmut_num/cover_num(ranges from 0 to 1).
An output example for the output pmat:
chr1 22525915 22525915 chr1_22525915_A. 968 0 0 0 A. A . 968 0 968 0.0
chr1 22525916 22525916 chr1_22525916_A. 968 0 0 0 A. A . 968 0 968 0.0
chr1 22525917 22525917 chr1_22525917_A. 974 0 0 0 A. A . 974 0 974 0.0
chr1 22525919 22525919 chr1_22525919_A. 977 0 0 0 A. A . 977 0 977 0.0
chr1 22525921 22525921 chr1_22525921_A. 975 0 0 0 A. A . 975 0 975 0.0
chr1 22525925 22525925 chr1_22525925_A. 982 0 0 0 A. A . 982 0 982 0.0
chr1 22525940 22525940 chr1_22525940_AG 897 128 0 0 AG A G 897 128 1025 0.12488
chr1 22525942 22525942 chr1_22525942_AG 968 59 0 0 AG A G 968 59 1027 0.05745
chr1 22525943 22525943 chr1_22525943_AG 1004 19 0 1 AG A G 1004 19 1023 0.01857
chr1 22525946 22525946 chr1_22525946_AG 996 27 0 0 AG A G 996 27 1023 0.02639
...
mpmat is one of the most critical intermediate files in the analysis workflow, merging tandem mutation sites within pmat into one signal region. The pmat-merge-ABE-v1.0.py script performs the conversion from pmat to mpmat. The following parameter settings are recommended.
# A-to-G on the Watson strand
python pmat-merge-ABE-v1.0.py \
-i mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20_A.pmat \
-f A -t G \
-r YOUR_PATH_TO/hg38_only_chromosome.fa \
-d 50 -D 100 --NoMutNumCutoff 1 --OmitTandemNumCutoff 1 \
-o mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.AG.mpmat \
2> mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.AG.mpmat.log &
# T-to-C on the Crick strand
python pmat-merge-ABE-v1.0.py \
-i mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20_T.pmat
-f T -t C \
-r YOUR_PATH_TO/hg38_only_chromosome.fa \
-d 50 -D 100 --NoMutNumCutoff 1 --OmitTandemNumCutoff 1 \
-o mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.TC.mpmat \
2> mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.TC.mpmat.log &Column explanation for mpmat:
chr_name: str, chromosome name of signal region, value like 'chr1', 'chr2' ...region_start: int, the start coordinate of the signal region, following a 1-based coordinate system.region_end: int, the end coordinate of the signal region, following a 1-based coordinate system.region_site_num: int, the total number of query nucleotides (e.g., A or T) contained within the region.region_contain_mut_site_num: int, the number of sites with query mutation (e.g., AG or TC, excluding SNPs) in the region.region_SNP_num: int, the number of sites with SNP in the region.region_site_index.list: str, a comma-separated string of site identifiers for each query site in the signal region. Each identifier is formatted as {chr_name}{site_coordinate}{query_base_type} (e.g., "chr1_10069_AG")to_base_num.list: str, a comma-separated string representing the number of reads supporting the query mutation at each corresponding site in <region_site_index.list>. The order of values matches the site order in <region_site_index.list> (e.g., "1,1" corresponds to 1 mutated read for each of two consecutive sites).cover_base_num.list: str, a comma-separated string representing the total number of reads covering each corresponding site in <region_site_index.list>. The order of coverage values is aligned with the site order in <region_site_index.list> (e.g., "12,13" means the first site has 12x coverage and the second has 13x coverage).mut_ratio.list: str, a comma-separated string representing the mutation ratio for each corresponding site in <region_site_index.list>.SNP_info.list: str, a comma-separated string indicating whether each site in <region_site_index.list> is a SNP (e.g., "False,False" for regions with no SNPs).tandem_info: str, a comma-separated '0' indicating tandem mutation site count in the region.
An output example for the output mpmat:
chr_name region_start region_end region_site_num region_mut_site_num region_SNP_mut_num region_site_index mut_base_num cover_base_num mut_ratio SNP_ann tandem_info
chr1 10063 10063 1 1 0 chr1_10063_AG 1 12 0.08333 False 0
chr1 10069 10074 2 2 0 chr1_10069_AG,chr1_10074_AG 1,1 12,13 0.08333,0.07692 False,False 0,0
chr1 10104 10104 1 1 0 chr1_10104_AG 1 16 0.0625 False 0
chr1 10566 10566 1 1 0 chr1_10566_AG 1 3 0.33333 False 0
chr1 10647 10647 1 1 0 chr1_10647_AG 1 2 0.5 False 0
chr1 10927 10927 1 1 0 chr1_10927_AG 1 2 0.5 False 0
chr1 13256 13256 1 1 0 chr1_13256_AG 1 27 0.03704 False 0
chr1 13491 13491 1 1 0 chr1_13491_AG 1 18 0.05556 False 0
chr1 13868 13868 1 1 0 chr1_13868_AG 3 7 0.42857 False 0
chr1 14542 14542 1 1 0 chr1_14542_AG 1 31 0.03226 False 0
...
The current version of mpmat contains all A-to-G and T-to-C mutation signal regions, requiring certain cutoffs in order to filter. This task is performed by the mpmat-select-ABE-v1.0.py command. The key parameters to focus on are --SiteMutNum (-m), --SiteCoverNum (-c), and --RegionPassNum. The cutoff values shown here are merely examples; typically, multiple combinations of cutoffs should be tested.
# A-to-G on the Watson strand
python mpmat-select-ABE-v1.0.py \
-i mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.AG.mpmat \
-o mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.AG.C4_M3_R2_T10.mpmat \
-f A -t G -m 3 -c 4 -r 0.01 --RegionPassNum 2 --RegionToleranceNum 10 --RegionMutNum 1 \
--InHeader True --OutHeader False &
# T-to-C on the Watson strand
python mpmat-select-ABE-v1.0.py \
-i mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.TC.mpmat \
-o mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.TC.C4_M3_R2_T10.mpmat \
-f T -t C -m 3 -c 4 -r 0.01 --RegionPassNum 2 --RegionToleranceNum 10 --RegionMutNum 1 \
--InHeader True --OutHeader False &Also signal regions in control samples MPMAT should be excluded:
# A-to-G on the Watson strand
bedtools sort \
-i mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.AG.C4_M3_R2_T10.mpmat \
-g YOUR_PATH_TO/hg38_only_chromosome.fa | bedtools intersect \
-b mpmat_select.AG/dIprofiling-ABE_control_rep1_hg38.MAPQ20.merge_d50_D100.AG.C4_M3_R2_T10.sort.mpmat -wa -v \
-sorted -g YOUR_PATH_TO/hg38_only_chromosome.fa.fai | uniq \
> mpmat_select.AG/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.AG.C4_M3_R2_T10.sort.Clean.mpmat &
# T-to-C on the Watson strand
bedtools sort \
-i mpileup/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.TC.C4_M3_R2_T10.mpmat \
-g YOUR_PATH_TO/hg38_only_chromosome.fa | bedtools intersect \
-b mpmat_select.TC/dIprofiling-ABE_control_rep1_hg38.MAPQ20.merge_d50_D100.TC.C4_M3_R2_T10.sort.mpmat -wa -v \
-sorted -g YOUR_PATH_TO/hg38_only_chromosome.fa.fai | uniq \
> mpmat_select.TC/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.TC.C4_M3_R2_T10.sort.Clean.mpmat &cat mpmat_select.merge/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.AG.C4_M3_R2_T10.Clean.sort.RmOverlap.mpmat \
mpmat_select.merge/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.TC.C4_M3_R2_T10.Clean.sort.RmOverlap.mpmat | \
grep -v chrM | bedtools sort -g YOUR_PATH_TO/hg38_only_chromosome.fa.fai | uniq \
> mpmat_select.merge/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.merge.C4_M3_R2_T10.Clean.sort.RmOverlap.RmChrM.sort.mpmat &Then find-significant-mpmat-ABE-v1.0.py command performs a comparison for those mutation signal regions between control sample and dI-profiling treatment sample by the Poisson statistical test. After this step, one can obtain tables containing information of signals regions (a mpmat file shown in the code part) and matched Poisson test results (a TSV table shown in the code part). The following are recommended parameter settings.
python find-significant-mpmat-ABE-v1.0.py \
-p 24 \
-i mpmat_select.merge/dIprofiling-ABE_treat_rep1_hg38.MAPQ20.merge_d50_D100.merge.C4_M3_R2_T10.Clean.sort.RmOverlap.RmChrM.sort.mpmat \
-o poisson_res/dIprofiling-ABE_treat_rep1__vs__ctrl.hg38.C4_M3_R2_T10.pvalue_table \
-c bam.bwa/dIprofiling-ABE_control_rep1_bwa_hg38_sort_rmdup.MAPQ20.bam \
-t bam.bwa/dIprofiling-ABE_treat_rep1_bwa_hg38_sort_rmdup.MAPQ20.bam \
-r YOUR_PATH_TO/hg38_only_chromosome.fa \
--query_mutation_type AG,TC \
--mpmat_filter_info_col_index -1 \
--mpmat_block_info_col_index -1 \
--region_block_mut_num_cutoff 2 \
--query_mut_min_cutoff 1 \
--query_mut_max_cutoff 8 \
--total_mut_max_cutoff 10 \
--other_mut_max_cutoff 6 \
--seq_reads_length 150 \
--lambda_method ctrl_max \
--poisson_method mutation 2> poisson_res/dIprofiling-ABE_treat_rep1__vs__ctrl.hg38.C4_M3_R2_T10.pvalue_table.log &Column explanation for pvalue_table:
chr_name: str, chromosome name of tested region, value like 'chr1', 'chr2' ...region_start: int, the start coordinate of the tested region, and the coordinate index is based on a 1-based scale.region_end: int, the end coordinate of the tested region, and the coordinate index is based on a 1-based scale.mpmat_index: str, a formatted string, which can be used as a key to index the whole table.region_site_num: int, number of sites (C or G) in the tested regionregion_block_site_num: int, number of sites that present a mutation signal in the control sample (A-to-G or T-to-C). The blocked sites are omitted in the enrichment test step.region_mut_site_num: int, number of sites with mutated signals in the treatment sample. Note, the blocked sites are not considered.region_site_index: str list, split by comma, list length is the same as <region_site_num>, and each item in this list is the site coordinate of the genome.region_block_state: str list, split by -, list length is the same as <region_site_num>, “B” means site is blocked, and "N" means site is not blocked.region_highest_site_index: str, coordinate of site with the highest dI-profiling signal.region_highest_site_mut_num: int, count of sequencing reads with tandem mutation info for the site with the highest dI-profiling signal.region_highest_site_cover_num: int, total count of sequencing reads for the site with the highest dI-profiling signal.region_highest_site_mut_ratio: float, range 0~1, mutation ratio, which equalsregion_highest_site_mut_num/region_highest_site_cover_numctrl_count: int, total count of sequencing reads in the control sample; if a read overlaps with the tested region, it will be counted.treat_count: int, total count of sequencing reads in treat sample.ctrl_mut_count: int, count of sequencing reads with A-to-G / T-to-C mutation info in control sample.treat_mut_count: int, count of sequencing reads with A-to-G / T-to-C mutation info in treat sample.ctrl_count.norm:, float, normalizedctrl_count, the default value equals count per million.treat_count.norm:, float, normalizedtreat_count, the default value equals count per million.ctrl_mut_count.norm:, float, normalizedctrl_mut_count, the default value equals count per million.treat_mut_count.norm:, float, normalizedtreat_mut_count, the default value equals count per million.count_info: meaning less in this version.log2_FC: float, log2 fold-change, which equals log2(treat_count.norm/ctrl_count.norm)log2_FC_mut: float, log2 fold-change, which equals log2(treat_mut_count.norm/ctrl_mut_count.norm)test_state: str, "TestOK" means the Poisson enrichment test works well.p_value: float, the p-value from the Poisson enrichment test.FDR: float, adjusted p-value with BH methods.
An output example for the enrichment significance test results pvalue_table:
chr_name region_start region_end mpmat_index region_site_num region_block_site_num region_mut_site_num region_site_index region_block_state region_highest_site_index region_highest_site_mut_num region_highest_site_cover_num region_highest_site_mut_ratio ctrl_count treat_count ctrl_mut_count treat_mut_count ctrl_count.norm treat_count.norm ctrl_mut_count.norm treat_mut_count.norm count_info log2_FC log2_FC_mut test_state p_value FDR
chr1 1467765 1467765 chr1_1467765_1467765 1 0 1 chr1_1467765_TC N chr1_1467765_TC 0 8 0.0 8 8 0 0 0.12133147645619574 0.09188187353939406 0.0 0.0 0,1 8,0 8,0 -0.4011016908785308 NA TestOK 0.6386278084328083 0.64141229048113
chr1 2197705 2197706 chr1_2197705_2197706 2 0 2 chr1_2197705_TC,chr1_2197706_TC N-N chr1_2197706_TC 6 17 0.35294117647058826 6 18 0 7 0.0909986073421468 0.20673421546363663 0.0 0.0803966393469698 0,1,2 6,0,0 11,3,4 1.1838608098426253 3.83987931615384 TestOK 0.015947120080862232 0.024188698946181787
chr1 2639405 2639406 chr1_2639405_2639406 2 0 2 chr1_2639405_AG,chr1_2639406_AG N-N chr1_2639405_AG 4 22 0.18181818181818182 2 23 0 4 0.030332869114048935 0.2641603864257579 0.0 0.04594093676969703 0,1,2 2,0,0 19,1,3 3.122460265178482 3.0325243940962356 TestOK 0.08280583178340417 0.10343565838688203
chr1 3476271 3476285 chr1_3476271_3476285 4 0 4 chr1_3476271_TC,chr1_3476276_TC,chr1_3476279_TC,chr1_3476285_TC N-N-N-N chr1_3476271_TC 12 24 0.5 6 28 0 16 0.0909986073421468 0.3215865573878792 0.0 0.18376374707878812 0,1,2,3,4 6,0,0,0,0 12,7,9,0,0 1.8212907304579171 5.032524394096235 TestOK 8.055007455896625e-05 0.00017839617739746512
chr1 3565744 3565749 chr1_3565744_3565749 2 0 2 chr1_3565744_AG,chr1_3565749_AG N-N chr1_3565744_AG 24 59 0.4067796610169492 4 67 0 27 0.06066573822809787 0.7695106908924253 0.0 0.31010132319545497 0,1,2 4,0,0 40,26,1 3.6649874995792415 5.787411896259704 TestOK 8.619981599539254e-08 2.659669536268094e-07
...
Select the significant dI-profiling signals according to the metrics in pvalue_table table with a certain criterion. In fact, dI-profiling across editors and biosamples yields heterogeneous patterns, each requiring a tailored threshold.
After finishing the selection step, one can obtain significantly enriched mutation signal regions of dI-profiling, indicating the off-target regions caused by adenine base editors. Run the mpmat-to-art-sgRNA-ABE-v1.0.py command if one want to further find sgRNA binding site for these regions. Although there are many parameters to set, the parameters provided here work in most cases.
python mpmat-to-art-sgRNA-ABE-v1.0.py \
-i final_list/293T-DI__ABE8e__ABEsite16.filter_candidates.IGVCheck.sort.mpmat \
-q GGGAATAAATCATAGAATCCNRG \
-r YOUR_PATH_TO/hg38_only_chromosome.fa \
-o final_list/293T-DI__ABE8e__ABEsite16.filter_candidates.IGVCheck.art \
-m region \
-e 40 \
--input_header False \
--mpmat_fwd_mut_type AG \
--mpmat_rev_mut_type TC \
--distal_index 1,10 \
--seed_index 11,20 \
--align_settings 5,-4,-24,-8 \
--PAM_type_penalty 0,4,8 \
--dna_bulge_penalty 24,8 \
--rna_bulge_penalty 24,8 \
--dna_bulge_cmp_weight 1,24 \
--rna_bulge_cmp_weight 1,24 \
--mismatch_cmp_weight 10,2 \
--dist_to_signal_penalty_k 0,0,0,0,0,0 \
--dist_to_signal_penalty_offset 12,0,0,0,0,12 > final_list/293T-DI__ABE8e__ABEsite16.filter_candidates.IGVCheck.align.log 2>&1 &Column explanation for art:
chrom: str, chromosome name of off-target region, value like 'chr1', 'chr2' ...start: int, the start coordinate of off-target region, and the coordinate index is based on a 1-based scale.end: int, the end coordinate of off-target region, and the coordinate index is based on a 1-based scale.region_index: str, the index of off-target region, used to uniquely identify different regions.align_chr_name: str, chromosome name of sgRNA aligned region, value like 'chr1', 'chr2' ...align_chr_start: int, the start coordinate of sgRNA aligned region, and the coordinate index is based on a 1-based scale.align_chr_end: int, the end coordinate of sgRNA aligned region, and the coordinate index is based on a 1-based scale.align_strand: str, the strand of sgRNA aligned region relative, value like '+' or '-'.align_dist_to_signal: int, the distance from the sgRNA aligned region to the off-target signal.PAM_type: str, the type of Protospacer Adjacent Motif (PAM) according to the sgRNA alignment, value like 'NGG', 'NAG' or 'NotNRG'.PAM_seq: str, the specific nucleotide sequence of the PAM, e.g., 'CGG' or 'TGG'.align_total_match: int, the total number of nucleotide that match between the off-target sequence and the on-target sequence in the alignment.align_total_mismatch: int, the total number of nucleotide that mismatch between the off-target sequence and the on-target sequence in the alignment.align_total_gap: int, the total number of gaps in the alignment between the off-target sequence and the on-target sequence.dna_bulge_num: int, the number of bulges (insertions) in the DNA sequence of the off-target sequence within the alignment relative to the on-target (sgRNA) sequence.rna_bulge_num: int, the number of bulges (insertions) in the sgRNA sequence within the alignment relative to the off-target sequence.head_mismatch: int, the number of mismatches in the head region (PAM-distal 10-bp) of the alignment between the off-target and on-target sequences.head_gap: int, the number of gaps in the head region (PAM-distal 10-bp) of the alignment between the off-target and on-target sequences.seed_mismatch: int, the number of mismatches in the seed region (PAM-proximal 10-bp) of the alignment between the off-target and on-target sequences.seed_gap: int, the number of gaps in the seed region (PAM-proximal 10-bp) of the alignment between the off-target and on-target sequences.align_penalty_score: int, the penalty score calculated for the alignment, considering mismatches and gaps (higher penalties indicate lower similarity).align_target_seq: str, the sequence of the off-target sgRNA alignment region, including gap characters (e.g., '-') to align with the on-target (sgRNA) sequence.align_info_state: str, the status of the alignment, value like '|' (match), '.' (mismatch) or '-' (gap).align_query_seq: str, the sequence of the on-target (sgRNA), including gap characters (e.g., '-') to align with the off-target sequence.
An output example for art:
chrom start end region_index align_chr_name align_chr_start align_chr_end align_strand align_dist_to_signal PAM_type PAM_seq align_total_match align_total_mismatch align_total_gap dna_bulge_num rna_bulge_num head_mismatch head_gap seed_mismatch seed_gap align_penalty_score align_target_seq align_info_state align_query_seq
...
chr5 82915258 82915271 chr5_82915258_82915271 chr5 82915249 82915271 - 0 NGG AGG 14 6 0 0 0 4 0 1 0 24 AGGAAAAGGACATAGAATCAAGG .||||.|...|||||||||...| GGGAATAAATCATAGAATCCNRG
chr5 128002824 128002826 chr5_128002824_128002826 chr5 128002815 128002837 - 0 NAG AAG 12 8 0 0 0 8 0 0 0 32 TCACTTTTTTCATAGAATCCAAG .....|...|||||||||||..| GGGAATAAATCATAGAATCCNRG
chr5 128125320 128125320 chr5_128125320_128125320 chr5 128125311 128125333 - 0 NGG TGG 13 7 0 0 0 7 0 0 0 28 GGTGTGGCCTCATAGAATCCTGG ||.......|||||||||||..| GGGAATAAATCATAGAATCCNRG
chr5 133205370 133205370 chr5_133205370_133205370 chr5 133205361 133205383 - 0 NGG AGG 15 5 0 0 0 4 0 1 0 20 AAGAAAAAGACATAGAATCCAGG ..|||.||..||||||||||..| GGGAATAAATCATAGAATCCNRG
chr5 136221877 136221893 chr5_136221877_136221893 chr5 136221868 136221890 - 0 NGG AGG 16 4 0 0 0 3 0 1 0 16 CCAAATAAACCATAGAATCCAGG ...||||||.||||||||||..| GGGAATAAATCATAGAATCCNRG
chr5 138461095 138461100 chr5_138461095_138461100 chr5 138461087 138461109 + 0 NGG AGG 12 8 0 0 0 7 0 1 0 32 TCTCTTGTAACATAGAATCCAGG .....|..|.||||||||||..| GGGAATAAATCATAGAATCCNRG
chr5 139403410 139403410 chr5_139403410_139403410 chr5 139403401 139403423 - 0 NGG TGG 15 5 0 0 0 4 0 1 0 20 GGGTAGGCATGATAGAATCCTGG |||.|...||.|||||||||..| GGGAATAAATCATAGAATCCNRG
chr5 143227663 143227670 chr5_143227663_143227670 chr5 143227658 143227679 + 0 NGG AGG 16 3 1 0 1 2 1 1 0 24 GGGAA-AGGCCATAGAATCCAGG |||||-|...||||||||||..| GGGAATAAATCATAGAATCCNRG
...
python plot-art-ABE-v1.0.py \
-i final_list/293T-DI__ABE8e__ABEsite16.final_list_subset.IGVCheck.art \
-o out_image/293T-DI__ABE8e__ABEsite16.final_list_subset.IGVCheck.align.pdf \
--align_seq GGGAATAAATCATAGAATCCTGG \
-k align_total_mismatch,seed_mismatch,align_total_gap,head_mismatch -r True,True,True,True \
-a align_coordinate,align_strand,align_total_mismatch,seed_mismatch,region_index &