v0.0.6

Author: yangao07

Makefile

@@ -11,14 +11,17 @@ endif
 
 # add -fno-tree-vectorize to avoid certain vectorization errors in O3 optimization
 # right now, we are using -O3 for the best performance, and no vectorization errors were found
-EXTRA_FLAGS = -Wall -Wno-unused-function -Wno-misleading-indentation -Wno-unused-variable -Wno-alloc-size-larger-than
+EXTRA_FLAGS = -Wall -Wno-misleading-indentation -Wno-unused-function #-Wno-unused-variable -Wno-alloc-size-larger-than
 
 # Define the version number
-LONGCALLD_VERSION =0.0.5
+VERSION=0.0.6
+# LONGCALLD_VERSION =0.0.6
 # Get the Git commit hash
 GIT_COMMIT := $(shell git rev-parse --short HEAD 2> /dev/null)
 ifneq ($(GIT_COMMIT),)
-	LONGCALLD_VERSION = 0.0.5-$(GIT_COMMIT)
+	LONGCALLD_VERSION = $(VERSION)-$(GIT_COMMIT)
+else
+	LONGCALLD_VERSION = $(VERSION)
 endif
 
 HTSLIB_DIR  = ./htslib
@@ -170,6 +173,6 @@ clean_all:
 clean_hts:
 	rm -f $(HTSLIB)
 clean_abpoa:
-	rm -f $(ABPOA_LIB)
+	rm -f $(ABPOA_LIB) $(ABPOA_DIR)/src/*.o
 clean_wfa2:
 	rm -f $(WFA2_LIB)

README.md

@@ -14,8 +14,9 @@
 ## Updates (pre-release v0.0.6)
 
 * Fix corrupted VCF output in v0.0.5
-* Low memory usage (especially when mosaic variant calling enabled)
-<!-- * Add [longdust](https://github.com/lh3/longdust) for long low-complexity regions -->
+* Fix missing MEI header in VCF output
+* Improved run time and memory usage (especially when mosaic variant calling enabled)
+* Add `--input-is-list` and `-X` to support multiple input BAM/CRAM files of the same sample for variant calling
 
 
 ## Getting Started
@@ -48,10 +49,12 @@ man ./longcallD.1
   - [Build from source](#build-from-source)
 - [Usage](#usage)
   - [Variant calling with PacBio HiFi/Nanopore long reads](#variant-calling-with-pacbio-hifinanopore-long-reads)
+  - [Variant calling with multiple input BAM/CRAM files of the same sample](#variant-calling-with-multiple-input-bamcram-files-of-the-same-sample)
   - [Low allele-frequency mosaic variant calling](#low-allele-frequency-mosaic-variant-calling)
   - [Region-specific variant calling](#region-specific-variant-calling)
-  - [Variant calling and output phased long reads](#variant-calling-and-output-phased-long-reads)
+  - [Variant calling and output phased (\& refined) long-read BAM/CRAM](#variant-calling-and-output-phased--refined-long-read-bamcram)
   - [Variant calling from remote files](#variant-calling-from-remote-files)
+- [Memory usage](#memory-usage)
 - [Acknowledgements](#acknowledgements)
 - [Contact](#contact)
 
@@ -107,14 +110,27 @@ longcallD call -t16 ref.fa hifi.bam > hifi.vcf         # default for PacBio HiFi
 longcallD call -t16 ref.fa ont.bam --ont > ont.vcf     # for ONT reads
 ```
 
+### Variant calling with multiple input BAM/CRAM files of the same sample
+You can provide multiple BAM/CRAM files of the same sample for variant calling using `--input-is-list` or `-X`:
+```
+longcallD call -t16 --input-is-list ref.fa bam_list.txt > sample.vcf
+# where bam_list.txt contains:
+# sample_part1.bam
+# sample_part2.bam
+# sample_part3.bam
+```
+or
+```
+longcallD call -t16 ref.fa sample_part1.bam -X sample_part2.bam -X sample_part3.bam > sample.vcf
+```
+
 ### Low allele-frequency mosaic variant calling
-With `-s`, longcallD will detect both germline and somatic/mosaic variants.
+With `-s`, longcallD will detect both germline and low-frequency somatic/mosaic variants.
 
 For each somatic/mosaic variant, a `SOMATIC` tag will be added to the INFO field in the output VCF.
 ```
 longcallD call -s -t16 ref.fa hifi.bam > hifi.vcf
 longcallD call -s -t16 ref.fa hifi.bam -T AluY_L1_SVA_cons_noPA.fa > hifi.vcf # add MEI information in INFO field
-longcallD call -s -t16 ref.fa ont.bam --ont > ont.vcf
 ```
 
 ### Region-specific variant calling
@@ -126,10 +142,10 @@ longcallD call -t16 ref.fa hifi.bam --region-file reg.bed > hifi_regs.vcf
 longcallD call -t16 ref.fa hifi.bam --autosome > hifi_autosome.vcf
 ```
 
-### Variant calling and output phased long reads
+### Variant calling and output phased (& refined) long-read BAM/CRAM
 ```
-longcallD call -t16 ref.fa hifi.bam --hifi -b hifi_phased.bam > hifi.vcf  # output phased HiFi reads (BAM tag: HP & PS)
-longcallD call -t16 ref.fa ont.bam --ont -b ont_phased.bam > ont.vcf      # output phased ONT reads (BAM tag: HP & PS)
+longcallD call -t16 ref.fa hifi.bam --hifi -b hifi_phased.bam > hifi.vcf                  # output phased HiFi reads (BAM tag: HP & PS)
+longcallD call -t16 ref.fa ont.bam --ont --refine-aln -b ont_phased_refined.bam > ont.vcf # output phased & refined ONT reads (BAM tag: HP & PS)
 ```
 ### Variant calling from remote files
 ```
@@ -138,6 +154,14 @@ bam=https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/HG002_NA2438
 longcallD call -t16 $ref $bam chr11:10,229,956-10,256,221 chr12:10,576,356-10,583,438 > hifi_regs.vcf
 ```
 
+## Memory usage
+As longcallD performs multiple-sequence alignment/re-alignment, which are memory-intensive, it usually uses more memory than other variant callers.
+The peak memory usage mainly depends on the number of threads (`-t/--threads`), the sequencing coverage, and the read length.
+For human genome sequencing data with ~40x coverage, longcallD typically uses around **1GB** (**HiFi**) or **2GB** (**ONT R10**) memory per thread for germline variant calling.
+
+If you encounter memory issues, you can use `--region-file` to limit the genomic regions being processed.
+Human genome region list excluding centromeres are provided [here](https://github.com/yangao07/longcallD/blob/main/anno/).
+
 ## Acknowledgements
 LongcallD is dependent on the following libraries, we are grateful to all the developers/maintainers:
 
@@ -146,7 +170,7 @@ LongcallD is dependent on the following libraries, we are grateful to all the de
 * [WFA](https://github.com/smarco/WFA2-lib): pairwise alignment
 * [edlib](https://github.com/Martinsos/edlib): fast sequence similarity calculation
 * [cgranges](https://github.com/lh3/cgranges): interval operations
-* [sdust](https://github.com/lh3/sdust) and [longdust](https://github.com/lh3/longdust): identify low-complexity regions
+* [sdust](https://github.com/lh3/sdust): identify low-complexity regions
 
 ## Contact
 

abPOA

@@ -0,0 +1,23 @@
+## Mobile element sequences
+File `AluY_L1_SVA_cons_noPA.fa` contains consensus sequences of three major
+human active mobile element families (AluY, L1HS and SVA) excluding the polyA
+tails. It can be used with the `-T` option to add mobile element insertion
+(MEI) information in the INFO field of somatic/mosaic variant calls.
+
+## Non-centromeric regions
+
+File `chm13v2.reg.nocen.bed` *excludes* the approximate locations of centromeric
+satellite repeats and acrocentric short arms in CHM13 v2.0, which was *manually*
+constructed based on the [official satellite][cen-sat] annotation, the
+[DNA-BRNN][dna-brnn] satellite annotation and the [minigraph pangenome
+graph][HPRC-mg] from the HPRC year-1 data.
+
+File `hs38.reg.nocen.bed` was constructed similarly from DNA-BRNN annotation and
+excludes regions where minigraph alignment faded.
+
+File `hs37.reg.nocen.bed` was constructed by running DNA-BRNN on GRCh37 (hs37d5.fa).
+
+[cen-sat]: https://s3-us-west-2.amazonaws.com/human-pangenomics/T2T/CHM13/assemblies/annotation/chm13v2.0_censat_v2.0.bed
+[dna-brnn]: https://github.com/lh3/dna-nn
+[HPRC-mg]: https://zenodo.org/records/10693675
+[zenodo]: https://zenodo.org/records/10963019

anno/README.md

@@ -0,0 +1,46 @@
+chr1	10000	121569169
+chr1	142292033	248377328
+chr2	10000	92250802
+chr2	94745067	242686752
+chr3	10000	90754701
+chr3	96465026	201095948
+chr4	10000	49655154
+chr4	55353192	193564945
+chr5	10000	46780042
+chr5	51012194	182035439
+chr6	10000	58236706
+chr6	61108390	172116628
+chr7	10000	60360644
+chr7	63764499	160557428
+chr8	10000	44193546
+chr8	46375080	146249331
+chr9	10000	44888599
+chr9	76744047	150607247
+chr10	10000	39583793
+chr10	41976237	134748134
+chr11	10000	50973358
+chr11	54526419	135117769
+chr12	10000	34543492
+chr12	37252490	133314548
+chr13	17558596	113556686
+chr14	12758411	101151492
+chr15	17744466	99743195
+chr16	10000	35784066
+chr16	52269756	96320374
+chr17	10000	23383372
+chr17	27621319	84266897
+chr18	10000	15591581
+chr18	21171235	80532538
+chr19	10000	24520766
+chr19	29819351	61697364
+chr20	10000	26333658
+chr20	33019590	66200255
+chr21	11356378	45080682
+chr22	15761065	51314926
+chrX	10000	2393410
+chrX	2395410	57769763
+chrX	60977195	153924834
+chrX	153926834	154249566
+chrY	10000	2457320
+chrY	2459320	62121809
+chrY	62123809	62450029