Reporting base counts for unedited sites

[eascarrunz]

Is it possible to report base counts at unedited sites with Reditools3?

Reditools2 did it in the non-strict mode, which seems gone in Reditools3. This information is useful for computing the Alu editing index.

[ahanden]

Yes, you can get total counts using Reditools3. Set the --min-edits option to zero and you'll get what you're looking for.

I'll make sure to add a note to our documentation about this too for other users.

[eascarrunz]

Thanks. I tried this command:

python -m reditools analyze ${bam} --reference ${genome} --strand ${strand_orientation} --output-file ${base_name}.site_edits_reditools3.txt --min-edits 0 --threads ${task.cpus} --verbose

But the output still skips sites without variants:

Region	Position	Reference	Strand	Coverage-q30	MeanQ	BaseCount[A,C,G,T]	AllSubs	Frequency	gCoverage-q30	gMeanQ	gBaseCount[A,C,G,T]	gAllSubs	gFrequency
21	6	C	-	1	37.00	[0, 0, 0, 1]	CT	1.00	-	-	-	-	-
21	510	A	-	5	37.00	[2, 0, 3, 0]	AG	0.60	-	-	-	-	-
21	741	A	-	10	37.00	[9, 0, 1, 0]	AG	0.10	-	-	-	-	-
21	1333	A	-	13	37.00	[3, 0, 10, 0]	AG	0.77	-	-	-	-	-
21	1466	C	-	9	37.00	[0, 0, 0, 9]	CT	1.00	-	-	-	-	-
21	1837	A	-	10	37.00	[9, 0, 1, 0]	AG	0.10	-	-	-	-	-

For comparison, the output from Reditools2:

Region	Position	Reference	Strand	Coverage-q30	MeanQ	BaseCount[A,C,G,T]	AllSubs	Frequency	gCoverage-q30	gMeanQ	gBaseCount[A,C,G,T]	gAllSubs	gFrequency
21	1	A	0	1	37.00	[1, 0, 0, 0]	-	0.00	-	-	-	-	-
21	2	C	0	1	37.00	[0, 1, 0, 0]	-	0.00	-	-	-	-	-
21	3	G	0	1	37.00	[0, 0, 1, 0]	-	0.00	-	-	-	-	-
21	4	G	0	1	37.00	[0, 0, 1, 0]	-	0.00	-	-	-	-	-
21	5	T	0	1	37.00	[0, 0, 0, 1]	-	0.00	-	-	-	-	-
21	6	C	0	1	37.00	[0, 0, 0, 1]	CT	1.00	-	-	-	-	-

[Xiezhouwei2024]

I also encountered the same problem. I set -me 0 as you said, but it couldn't be solved. If the non-editing sites are deleted, REDInet will not be used for subsequent predictions. Currently it can only be run through reditools2, but the 2nd generation software is slower.

[ahanden]

Can you confirm which release/branch of REDI3 you're using?

[Xiezhouwei2024]

I think it's the latest version～
sh-4.2$ pip list
Package Version

---

pip 25.2
pysam 0.23.3
REDItools3 3.4
setuptools 80.9.0
sortedcontainers 2.4.0
wheel 0.45.1

[ahanden]

For v3.4, you also need to set --variants to "all". If you can post the full function call you're using, I may be able to see if anything else is amiss. Otherwise, I'm hoping to get v3.5 out by the end of the weekend which I have been able to get non-edit sites without issue.

[Xiezhouwei2024]

Thank you very much for your prompt reply, my code is as follows:

python -m reditools analyze \
/data/input/Files/ResultData/Workflow/W202505080000130/02.mergeBam/SC_T21E5_Cerebellum_1_NB.filter.sorted.bam \
-r /data/input/Files/chenjin/W202504210003221/reference/WGS_hg38.fa \
-o /data/work/test_new.txt \
-s 0 -C -V -t 10 -men 5 -me 0 -v all -q 25 -bq 25 

1. By the way,I don't fully understand the meaning of some of these parameters, and I hope you can provide a more detailed explanation of the parameters in the new version.

2. Taking my code as an example, I found that -men, -qand -bq have an impact on the results. I set -men 5 based on the suggestions of other users, and I did not fully understand its significance.

3. At the same time, I found that increasing -t did not speed up the running of the software. Why is this? Because the running time of retaining non-edited sites is relatively long, I am currently using reditools2 to run, and one chromosome needs to be run for an afternoon. Is there any way to speed up the process in reditools3?

[ahanden]

First thing, I was able to get the next major update out yesterday. You can install reditools v3.5 with this: pip install reditools --upgrade.

The issue with your current call is the -men 5. This option works differently in REDI3 vs REDI2, and does need to be addressed. For REDI3, the -men option identifies bases where there are multiple alts. Dropping that parameter will get you all variants. I am working on a detailed documentation as well as a tutorial that should be done by the end of December.

The -q option filters by MAPQ, which is the alignment score produced by whatever aligner you chose to use.

The -bq option filters by PHRED score, which comes from your FASTQ files. Depending on whether or how you trimmed your FASTQ files before alignment, this option may not be necessary for your analysis.

The -t option does spawn multiple threads, but it does have an effective cap. Each thread will receive one contig of your alignment file. If you are working with a genome of 24 chromosomes, you will not get any speedup beyond -t 24.

The way to maximize the speedup of reditools is to make sure your BAM file has MD tags. If you give reditools a bam file with MD tags and no reference fasta file, you should see a speedup of at least 20x.

[Xiezhouwei2024]

Thank you vvvvvery much, I ran it successfully using V3.5, the code is as follows:

#Reditools3
python -m reditools analyze \
/data/users/liyuejiao/online/project/t21/rna_editing/dnbc4tools/T21E6_Midbrain/T21E6_Midbrain/outs/anno_decon_sorted.bam \
-r /data/work/REDInet_all/REDInet/Package/Utilities/GRCh38.primary_assembly.genome.fa \
-o /data/work/REDItools3-3.5-stable/T21_result/T21E6_Midbrain.txt \
-s 0 -C -V -t 24 -q 25 -me 0  -g chr21

1. I deleted -men 5 in reditools3 according to your suggestion. My understanding is that it will filter the bases of less than 5 variants, resulting in the inability to output complete results. Runs successfully and gets results containing non-edited sites.

After comparing 2 and 3, I finally decided to use Reditools2 because it is easier to run. I still retain -men 5 in reditools2. Is this correct?
As you said, this is human short-read data. -t 24 is the maximum value in Reditools3. The result of one chromosome is about 3-6G (the running time is about 2-6h). If the entire sample of 24 chromosomes is run at the same time, the speed is still very slow.

My current strategy is to set -g chr1, 2, 3... to run a single chromosome, with one chromosome as one task, split into 24 tasks in total, and finally run REDInet to get the RNA editing results of each chromosome, finally merge them. My final code is as follows：

#Reditools2
python3 /data/work/Reditools2_all/reditools2.0-master/src/cineca/reditools.py \
-f /data/users/liyuejiao/online/project/t21/rna_editing/dnbc4tools/T21E6_Midbrain/T21E6_Midbrain/outs/anno_decon_sorted.bam \
-r /data/input/Files/chenjin/W202504210003221/reference/WGS_hg38.fa \
-o /data/work/Reditools2_all/Reditools2_result/T21E6_chr1.txt \
-l 1 -q 25 -bq 25 -men 5 -C -g chr1

[ahanden]

The -l 1 option ought to exclude any positions where you have no edit events. If that is what you intended, then that all looks great.

Depending on how much computing power you have, you could split up your analysis further. This is a version of a script I often used with REDI2 using a super computing cluster (modified to use your parameters)

#!/bin/bash
window_size=10000000
bam="/data/users/liyuejiao/online/project/t21/rna_editing/dnbc4tools/T21E6_Midbrain/T21E6_Midbrain/outs/anno_decon_sorted.bam"
stem=`basename $bam | sed 's/\.bam//'`
let idx=0
samtools view -H $bam | \
grep '^@SQ' | \
while read seq_line; do
        contig=`echo $seq_line | awk '{print $2}' | sed 's/SN://'`
        length=`echo $seq_line | awk '{print $3}' | sed 's/LN://'`
        for start in `seq 0 $window_size $length`; do
                ((idx++))
                idx_str=`printf %03d $idx`
                output="${stem}_${idx_str}.txt"
                stop=`expr $start + $window_size`
                window="$contig:$start-$stop"
                python3 /data/work/Reditools2_all/reditools2.0-master/src/cineca/reditools.py \
                        -f $bam \
                        -r /data/input/Files/chenjin/W202504210003221/reference/WGS_hg38.fa \
                        -o /data/work/Reditools2_all/Reditools2_result/$output \
                        -l 1 -q 25 -bq 25 -men 5 -C -g $window
        done
done

[Xiezhouwei2024]

Regarding '-l 1' what I want to retain is the information of all sites, including non-edited sites, because subsequent REDInet runs require continuous nucleotide information. So I guess I need to remove '-l 1', but the default value of -l is 1

sh-4.2$ grep "min-column-length" /data/work/Reditools2_all/reditools2.0-master/src/cineca/reditools.py
    parser.add_argument('-l', '--min-column-length', type=int, default=1, help='The minimum length of editing column (per position). Positions whose columns have length below this value will not be included in the analysis.')

When I use '-l 1' the result is that I can still run REDInet normally. Here is my running report:

Final Report:
Inference on T21E6_chr2.txt start.
Inference on T21E6_chr2.txt end. Elapsed time 0:00:05.654936.
[2025-11-24 14:29:49.977580] T21E6_chr2.txt processing start.
Total evaluated rows in T21E6_chr2.txt: 118220947
Total discarded sites in T21E6_chr2.txt: 14624.
Total candidates sites in T21E6_chr2.txt: 297968.
Candidates sites identification in T21E6_chr2.txt end. Elapsed time: 0:03:30.421747.
Total sequence with eccessive missing nucleotides in T21E6_chr2.txt: 296167.
Total extracted sites in T21E6_chr2.txt: 1801.
Features extraction in T21E6_chr2.txt end. Elapsed time 0:25:39.762869

You can see 297968 candidates sites but 296167 is missing nucleotides. Is it because of ‘-l 1‘?