Init demuxalot

[mschecht]

Add demuxalot demultiplexing to workflow

This issue tracks the addition of demuxalot to the workflow. This will be the first extension of the snakemake pipeline beyond Cell Ranger into downstream single-cell preprocessing steps - expect some bumps along the way!

Background

Preprint: Rogozhnikov et al., 2021

Implementation Notes

1. Two-stage approach see issue #21]:

   • Start with boilerplate demultiplexing (use existing SNPs)
   • Later: add SNP identification workflow

2. Test strategy:

   • Implement boilerplate workflow first
   • Validate on test dataset before expanding

Helper Script: Parse demuxalot Output

def process_demuxalot_results(base_dir, confidence_threshold=0.9):
    """
    Parse demuxalot posterior probabilities into assignment dataframe.
    
    Args:
        base_dir: Directory containing demuxalot output subdirectories
        confidence_threshold: Minimum probability for confident assignment (default: 0.9)
    
    Returns:
        DataFrame with columns: BARCODE, assignment, prob_max, assignment_status_demuxalot, sample
    """
    all_results = []
    
    for subdir in sorted(os.listdir(base_dir)):
        subdir_path = os.path.join(base_dir, subdir)
        if not os.path.isdir(subdir_path):
            continue
        
        prob_path = os.path.join(subdir_path, "posterior_probabilities.tsv.gz")
        if not os.path.exists(prob_path):
            continue
            
        # Load posterior probabilities
        prob_df = pd.read_csv(prob_path, sep="\t").set_index("BARCODE")
        
        # Assign max genotype and compute confidence
        assignments = prob_df.idxmax(axis=1)
        max_probs = prob_df.max(axis=1)
        
        results_df = pd.DataFrame({
            "BARCODE": prob_df.index,
            "assignment": assignments,
            "prob_max": max_probs,
        })
        
        # Filter by confidence threshold
        results_df.loc[results_df["prob_max"] < confidence_threshold, "assignment"] = "unassigned"
        
        # Classify assignment type
        def classify(x):
            if x == "unassigned":
                return "unassigned"
            elif "+" in x.lower() or x == "doublet":
                return "doublet"
            else:
                return x
        
        results_df["assignment_status_demuxalot"] = results_df["assignment"].apply(classify)
        
        # Add sample identifier
        results_df["sample"] = subdir.replace("_doublet_prior", "")
        all_results.append(results_df)
    
    # Concatenate all results
    final_df = pd.concat(all_results, ignore_index=True)
    return final_df

# Example usage
demuxalot_results_df = process_demuxalot_results(
    os.path.join(WD, "04_DEMUX_genotypes_lists_demuxalot")
)

TODO

• Add demuxalot to environment.yml
• Expand test dataset workflow to this step!
  • Need a BAM file from cellranger to compute genotype VCF files
• Include demuxalot commands in config generator
• Create boilerplate demultiplexing rule
• Integrate parsing script into workflow
• Add SNP identification workflow (future)
• Documentation

[mschecht]

First step, installation as seen on their README:

pip install demuxalot

Then export new conda yamls!

[mschecht]

This is how we are making test VCFs from the cellranger test data so we can develope demux tools.

Step 1. Run test dataset but make intermediate bam files

Modify the test_config_gex.yaml

$ cat 00_TEST_DATA_GEX/test_config_gex.yaml
cat 00_TEST_DATA_GEX/test_config_gex.yaml
project_name: test_gex
output_dir: test_output_gex
resources:
  mem_gb: 64
  tmpdir: ''
directories_suffix: none
cellranger_gex:
  enabled: true
  reference: /project/lbarreiro/SHARED/PROGRAMS/cellranger-9.0.1/external/cellranger_tiny_ref
  libraries: 00_TEST_DATA_GEX/libraries_list_gex.tsv
  chemistry: auto
  normalize: none
  create-bam: true # TURN ON
  threads: 10
  mem_gb: 64
doublet_detection:
  enabled: false
  method: scrublet
  scrublet:
    expected_doublet_rate: 0.06
    min_counts: 2
    min_cells: 3
celltype_annotation:
  enabled: false
  method: celltypist
  celltypist:
    model: Immune_All_Low.pkl
    majority_voting: false

Run the test data set to generate the BAM files.

snakemake-run-cellranger run --config-file 00_TEST_DATA_GEX/test_config_gex.yaml --cores 1

Make consensus VCF with bcftools - @danielsarj figured this out!

cd cellranger-snakemake/tests

echo "tests/test_output_gex/01_CELLRANGERGEX_COUNT/1_L001/outs/possorted_genome_bam.bam" > bamfiles.txt
echo "tests/test_output_gex/01_CELLRANGERGEX_COUNT/1_L002/outs/possorted_genome_bam.bam" >> bamfiles.txt

FASTA="/path/to/cellranger-9.0.1/external/cellranger_tiny_ref/fasta/genome.fa"

bcftools mpileup -Ou -f $FASTA --bam-list bamfiles.txt | bcftools call -mv -Oz -o L001_L002.vcf.gz 
tabix -p vcf L001_L002.vcf.gz

Test demux tools, here's an example with cellsnp-lite + demuxalot

# Demux with cellsnp-lite + vireo
cellsnp-lite \
    -s test_output_gex/01_CELLRANGERGEX_COUNT/1_L001/outs/possorted_genome_bam.bam \
    -b test_output_gex/01_CELLRANGERGEX_COUNT/1_L001/outs/filtered_feature_bc_matrix/barcodes.tsv.gz \
    -O cellsnp_output_L001 \
    -R tests/L001_L002.vcf.gz \
    -p 4 \
    --minMAF 0.0 \
    --minCOUNT 1 \
    --UMItag Auto \
    --cellTAG CB \
    --gzip

pip install -U vireoSNP
vireo -c cellsnp_output_L001 -N 2 -o vireo_output_L001_k2

$ ls vireo_output_L001_k2
donor_ids.tsv  fig_GT_distance_estimated.pdf  GT_donors.vireo.vcf.gz  _log.txt  prob_doublet.tsv.gz  prob_singlet.tsv.gz  summary.tsv

Notes:

• Need to decide if we are going to version on git or put the test vcf on zenodo