• CNVeil
  • Install CNVeil
  • Run CNVeil
    • Step 1: Extract read depth matrix
    • Step 2: Infer total copy number
    • Step 3: Infer allele-specific copy number
    • Step 4: Infer haplotype-specific copy number
• Benchmark
  • Prepare scDNA-seq data
  • Run Tools

You can install CNVeil through github

git clone git@github.com:maiziezhoulab/CNVeil.git

To set up the environment, you should have conda installed on your machine. Then simply run

# create conda env CNVeil_py3
conda env create -f CNVeil/env_yamls/environment_py3.yml 
# create conda env CNVeil_py2
conda env create -f CNVeil/env_yamls/environment_py2.yml
# create conda env CNVeil_R451
conda env create -f CNVeil/env_yamls/environment_R451.yml

You will have 3 virtual environments installed (CNVeil_py3, CNVeil_py2 and CNVeil_R451). And they will be used in diferent steps. You will need to load the corresponding virtual environment before you run each step in CNVeil pipeline.

This step bins reads from the input BAM files, performs GC correction, and reformats the corrected bin-level read depth into a matrix for downstream CNVeil analysis.

Script:

bin/CNVeil_extract_read_depth.py

Required conda environment:

conda activate CNVeil_R451

For real data:

python bin/CNVeil_extract_read_depth.py \
  --merged_bam_file /path/to/merged.bam \
  --split_bam_dir /path/to/split_bams \
  --work_dir /path/to/work_dir \
  --reftype hg38 \
  --num_threads 20

For simulation data with a cell-node mapping file:

python bin/CNVeil_extract_read_depth.py \
  --merged_bam_file /path/to/merged.bam \
  --split_bam_dir /path/to/split_bams \
  --work_dir /path/to/work_dir \
  --reftype hg38 \
  --num_threads 20 \
  --cell_node /path/to/cell_node.tsv

The script creates the following output directory:

<work_dir>/Total_CN/

Main output:

<work_dir>/Total_CN/RDmatrix.csv

This is the final read-depth matrix used by downstream CNVeil steps.

This step takes the read-depth matrix generated from Step 1, and outputs a total copy-number matrix for downstream CNVeil analysis.

Script:

bin/CNVeil_total_CN.py

Required conda environment:

conda activate CNVeil_py3

This script expects the Step 1 output directory structure:

<work_dir>/Total_CN/

Required input file:

<work_dir>/Total_CN/RDmatrix.csv

For real 10x data:

python bin/CNVeil_total_CN.py \
  --work_dir /path/to/work_dir \
  --reftype hg38 \
  --seq_type <paired-end or single-end> \
  --data_type 10x \
  --sample SAMPLE_NAME \
  --num_threads 20

For other single-cell sequencing data (non-10x):

python bin/CNVeil_total_CN.py \
  --work_dir /path/to/work_dir \
  --reftype hg38 \
  --seq_type <paired-end or single-end> \
  --data_type other \
  --sample SAMPLE_NAME

For simulation data with a cell-node mapping file:

python bin/CNVeil_total_CN.py \
  --work_dir /path/to/work_dir \
  --reftype hg38 \
  --seq_type <paired-end or single-end> \
  --data_type 10x \
  --sample SAMPLE_NAME \
  --cell_node /path/to/cell_node.tsv

All outputs are written under:

<work_dir>/Total_CN/

Main total copy-number matrix:

<work_dir>/Total_CN/CNV_<sample>.csv

This file has cells as rows and genomic bins as columns.

Transposed copy-number matrix:

<work_dir>/Total_CN/CNV_<sample>.csv.T

This file has genomic bins as rows and cells as columns.

This step uses the total copy-number results from Step 2 together with SNP allele-count information to infer allele-specific copy number (ASCN). You need to first extract per cell allele count before running the main script.

To get allele-count information, follow the instruction here.

After that, you will get var_all.rds in the output folder and you need to process it using bin/load_snp.R.

Activate the R environment before running:

conda activate CNVeil_R451

Run:

Rscript load_snp.R <rds_dir> <work_dir> <sample> <num_threads>

Example:

Rscript load_snp.R \
  /path/to/rds_dir \
  /path/to/work_dir \
  TN1 \
  8

Required input files in <rds_dir>:

var_all.rds
ref_all.rds
alt_all.rds

Required total CN file:

<work_dir>/Total_CN/CNV_<sample>.csv

This step outputs files to:

<work_dir>/Allele_CN/allele_count_by_cell/

Main outputs:

<sample>_baf.tsv
<sample>_baf_chr1.tsv
...
<sample>_baf_chr22.tsv
<sample>_pos_baf_chr1.tsv.gz
...
<sample>_pos_baf_chr22.tsv.gz

The *_baf_chr*.tsv files contain cell-level allele counts:

chrom    pos    cell    ref    alt

The *_pos_baf_chr*.tsv.gz files contain position-level BAF features:

chrom    pos    ref_sum    alt_sum    n_cells    n_obs    depth    baf    mbaf

Normal cells are inferred from average total CN:

|avg_CN - 2| <= 0.1

Normal cells are excluded when computing position-level BAF files.

Afterwards, you may run bin/CNVeil_Allele_CN.py to generate allele-specifc CN.

Script:

bin/CNVeil_allele_CN.py

Required conda environment:

conda activate CNVeil_py3

This step expects the CNVeil working directory to already contain the total CN result from Step 2:

<work_dir>/Total_CN/CNV_<sample>.csv

It also requires a directory containing SNP/RDS-derived input files:

<var_dir>/

The SNP loader expects this directory to contain:

var_all.rds
ref_all.rds
alt_all.rds

Basic usage:

python bin/CNVeil_allele_CN.py \
  --work_dir /path/to/work_dir \
  --var_dir /path/to/var_rds_dir \
  --ref_type hg38 \
  --num_threads 22

With a cell-name converter:

python bin/CNVeil_allele_CN.py \
  --work_dir /path/to/work_dir \
  --var_dir /path/to/var_rds_dir \
  --ref_type hg38 \
  --convert /path/to/collected_barcodes_SAMPLE.txt \
  --num_threads 22

When running for 10x data:

python bin/CNVeil_allele_CN.py \
  --work_dir /path/to/work_dir \
  --var_dir /path/to/var_rds_dir \
  --ref_type hg38 \
  --convert /path/to/collected_barcodes_SAMPLE.txt \
  --dtype 10x \
  --num_threads 22

The script writes allele-specific CN results under:

<work_dir>/Allele_CN/

Important outputs include:

<work_dir>/Allele_CN/ascn_seg.csv

Allele-specific copy number inferred at the cell/segment level.

This step uses the allele-specific copy-number results from Step 3 and phases allele-specific states along the genome to infer haplotype-specific copy number (HAP-CN).

Script:

bin/CNVeil_hap_CN.py

Required conda environment:

conda activate CNVeil_py2

This step expects the CNVeil working directory to already contain the allele-specific CN result from Step 3:

<work_dir>/Allele_CN/ascn_seg.csv

This file contains allele-specific copy number states per cell and segment, formatted as:

chr,start,end,cell1,cell2,...
chr1,100000,200000,2|1,1|2,...

Each entry A|B represents the allele-specific copy numbers for a given segment and cell.

Basic usage:

python2.7 bin/CNVeil_hap_CN.py \
  --ascn_file /path/to/work_dir/Allele_CN/ascn_seg.csv \
  --out_dir /path/to/work_dir/Haplotype_CN 

The script writes haplotype-specific CN results under:

<work_dir>/Haplotype_CN/

Important outputs include:

Haplotype-specific copy number inferred at the cell/segment level with consistent allele orientation across the genome.

<work_dir>/Haplotype_CN/phasing_output.pkl

Serialized phased copy-number object for downstream analysis.

Here is an example of generating the bam files from fastq files downloaded from NCBI cell by cell. The patient is named as 'KTN126'.

patient="KTN126"
file=YOUR_PATH_FOR_SRR_LIST
srapath=YOUR_PATH_TO_STORE_SRA_FILE
fastqpath=YOUR_PATH_TO_STORE_FASTQ_FILE
mkdir ${fastqpath}

 while read line; do
     echo $line
     YOUR_PATH_TO_SRATOOLKIT/bin/fastq-dump -I --split-files $srapath$line"/"$line".sra" 
     -O $fastqpath$line
 done < $file

Generate the fastq files from the. sra file with sra toolkit.

ml Intel/2017.4.196 BWA/0.7.17
ml GCC/11.3.0 SAMtools/1.18
sampath=${srapath}/sampath_SraAccList${patient}/
align_dir=${srapath}/sampath_SraAccList${patient}/
readgroup_path=${srapath}/bampath_SraAccList${patient}_readgroup/
dedup_path=${srapath}/bampath_SraAccList${patient}_dedup/

mkdir ${sampath} ${readgroup_path} ${dedup_path}
 
while read line; do
    echo $line
    cd $align_dir 
     bwa mem -M -t 8 YOUR_PATH_TO_REF_GENOME_FILE $fastqpath/$line/"$line"_1.fastq $fastqpath/$line/"$line"_2.fastq > $align_dir/"$line".sam
     samtools view -bS $align_dir/"$line".sam > $align_dir/"$line".bam 
     java -jar YOUR_PATH_TO_PICARD/picard.jar SortSam -I $align_dir/"$line".bam -SORT_ORDER coordinate -O $align_dir/"$line"_hg38.sorted.bam 
     java -jar YOUR_PATH_TO_PICARD/picard.jar AddOrReplaceReadGroups -I $align_dir/"$line"_hg38.sorted.bam -O $readgroup_path/"$line"_hg38.sorted.rg.bam -RGID $readgroup_path/"$line"_hg38/ -RGLB NAVIN_Et_Al -RGPL ILLUMINA -RGPU machine -RGSM $readgroup_path/"$line"_hg38/
     java -jar YOUR_PATH_TO_PICARD/picard.jar MarkDuplicates -REMOVE_DUPLICATES true -I $readgroup_path/"$line"_hg38.sorted.rg.bam -O $dedup_path/"$line"_hg38.sorted.rg.dedup.bam -METRICS_FILE $dedup_path/"$line"_hg38.sorted.rg.dedup.metrics.txt -PROGRAM_RECORD_ID MarkDuplicates -PROGRAM_GROUP_VERSION null -PROGRAM_GROUP_NAME MarkDuplicates
     java -jar YOUR_PATH_TO_PICARD/picard.jar BuildBamIndex -I $dedup_path/"$line"_hg38.sorted.rg.dedup.bam -O $dedup_path/"$line"_hg38.sorted.rg.dedup.bai
done < $file

Generate the bam files after alignment.

Next, we generate the bed files with hg19/hg38 reference.

ml GCC/8.2.0 BEDTools/2.28.0
bedpath_hg38=${srapath}/bedpath_SraAccList${patient}_hg38
bedpath_hg19=${srapath}/bedpath_SraAccList${patient}_hg19
mkdir ${bedpath_hg38} ${bedpath_hg19}
while read line; do
    echo $line
    bamToBed -i $dedup_path/$line'_hg38.sorted.rg.dedup.bam' > \
        ${bedpath_hg38}/$line'_grch38_rmdup.bed'
done < $file

while read line; do
    echo $line
    YOUR_PATH_TO_LIFTOVER/liftOver \
        ${bedpath_hg38}/$line'_grch38_rmdup.bed' \
        'YOUR_PATH_TO_LIFTOVER_CHAINFILE/liftOver_ChainFiles/hg38ToHg19.over.chain' \
        ${bedpath_hg19}/$line'_hg19_rmdup.bed' \
        ${bedpath_hg19}/$line'_hg19_unlifted_rmdup.bed'
done < $file

Zipped the bed file, which is recommended for the Ginkgo web platform to analyze.

cd ${bedpath_hg19}