2 Cell Ranger and Demultiplexing
Each multiplexed sample contains cells from two donors (one PASC, one Healthy) pooled together before library preparation. The processing pipeline has three stages: (1) Cell Ranger ARC jointly processes the paired RNA and ATAC libraries to produce a common feature-barcode matrix; (2) SNP-based demultiplexing assigns each cell barcode to a donor; (3) HLA-based approaches independently confirm which donor carries PASC versus Healthy status.
2.1 Cell Ranger ARC
Raw FASTQ files were processed with Cell Ranger ARC v2.0.2, which is designed for
10x Multiome (joint RNA + ATAC) libraries. Each sample has four ATAC FASTQ sets and one
Gene Expression FASTQ set, described in a libraries CSV file that Cell Ranger ARC uses
to link the two modalities. The --reference flag points to a custom GRCh38 multiome
reference. Outputs include a joint filtered feature-barcode matrix (used for all
downstream R analysis), a fragment file for ATAC peak calling, and sorted BAM files for
both the RNA and ATAC libraries.
cat << EOF >libraries_B.csv
fastqs,sample,library_type
rawdata,ATAC-Sample-B-1,Chromatin Accessibility
rawdata,ATAC-Sample-B-2,Chromatin Accessibility
rawdata,ATAC-Sample-B-3,Chromatin Accessibility
rawdata,ATAC-Sample-B-4,Chromatin Accessibility
rawdata,cDNA-Sample-B-1,Gene Expression
EOF
~/Programs/cellranger-arc-2.0.2/cellranger-arc count --id=cellranger_SC_ATAC_covid_B \
--reference=GRCh38_ScATAC \
--libraries=libraries_B.csv \
--localcores=16
cat << EOF >libraries_C.csv
fastqs,sample,library_type
rawdata,ATAC-Sample-C-1,Chromatin Accessibility
rawdata,ATAC-Sample-C-2,Chromatin Accessibility
rawdata,ATAC-Sample-C-3,Chromatin Accessibility
rawdata,ATAC-Sample-C-4,Chromatin Accessibility
rawdata,cDNA-Sample-C-1,Gene Expression
EOF
~/Programs/cellranger-arc-2.0.2/cellranger-arc count --id=cellranger_SC_ATAC_covid_C \
--reference=GRCh38_ScATAC \
--libraries=libraries_C.csv \
--localcores=16
cat << EOF >libraries_D.csv
fastqs,sample,library_type
rawdata,ATAC-Sample-D-1,Chromatin Accessibility
rawdata,ATAC-Sample-D-2,Chromatin Accessibility
rawdata,ATAC-Sample-D-3,Chromatin Accessibility
rawdata,ATAC-Sample-D-4,Chromatin Accessibility
rawdata,cDNA-Sample-D-1,Gene Expression
EOF
~/Programs/cellranger-arc-2.0.2/cellranger-arc count --id=cellranger_SC_ATAC_covid_D \
--reference=GRCh38_ScATAC \
--libraries=libraries_D.csv \
--localcores=162.2 Demultiplexing using SNPs - identify Donor
Because each sample is a pool of two donors, cell barcodes must be assigned to individual donors before any biological analysis. SNP-based demultiplexing exploits naturally occurring germline genetic variants present in the RNA reads: cells from different donors carry different alleles at heterozygous SNP positions, so clustering cells by their allele profiles separates the two donors without any prior genotype information.
scSNPdemux implements this approach using
cellsnp-lite to genotype each cell at a panel of common variants
(genome1K.phase3.SNP_AF5e2 — ~5% minor allele frequency SNPs from the 1000 Genomes
Phase 3 reference panel, hg38), then runs Vireo to assign cells to donors. The -n 2
flag specifies that two donors are present in each pool. Cells that cannot be confidently
assigned are labelled unassigned; cells whose allele profiles are inconsistent with
either donor singlet are labelled doublet. Both categories are excluded from all
downstream analysis.
# https://github.com/wkljohn/scSNPdemux
git clone https://github.com/wkljohn/scSNPdemux.git
# install and create conda environment using instructions from the GitHub README
wget https://sourceforge.net/projects/cellsnp/files/SNPlist/genome1K.phase3.SNP_AF5e2.chr1toX.hg38.vcf.gz
conda activate demuxEnvironment
./demux.sh -b ../cellranger_SC_ATAC_covid_B/outs/gex_possorted_bam.bam \
-i ../cellranger_SC_ATAC_covid_B/outs/filtered_feature_bc_matrix/ \
-p genome1K.phase3.SNP_AF5e2.chr1toX.hg38.vcf.gz \
-o cellranger_SC_ATAC_covid_B.demux -n 2
./demux.sh -b ../cellranger_SC_ATAC_covid_C/outs/gex_possorted_bam.bam \
-i ../cellranger_SC_ATAC_covid_C/outs/filtered_feature_bc_matrix/ \
-p genome1K.phase3.SNP_AF5e2.chr1toX.hg38.vcf.gz \
-o cellranger_SC_ATAC_covid_C.demux -n 2
./demux.sh -b ../cellranger_SC_ATAC_covid_D/outs/gex_possorted_bam.bam \
-i ../cellranger_SC_ATAC_covid_D/outs/filtered_feature_bc_matrix/ \
-p genome1K.phase3.SNP_AF5e2.chr1toX.hg38.vcf.gz \
-o cellranger_SC_ATAC_covid_D.demux -n 2The output donor_ids.tsv file from each run maps every cell barcode to donor0,
donor1, doublet, or unassigned. These labels are anonymous — they do not carry
biological meaning on their own. The PASC/Healthy identity of each donor number is
resolved in the HLA-based steps below.
2.3 Demultiplexing Using HLA Haplotypes - identify PASC or Healthy
SNP demultiplexing separates the two donors but assigns them arbitrary labels (donor0,
donor1). To determine which donor is PASC and which is Healthy, HLA allele expression
is used as an orthogonal identifier. Each individual carries a unique combination of HLA
alleles — alleles that are private to one donor can be used as molecular barcodes to
confirm identity and link the anonymous donor labels to the known clinical metadata.
A custom Cell Ranger reference is built from FASTA and GTF files describing the
donor-specific HLA allele sequences. Cell Ranger is then run in standard count mode
against this custom reference, treating each HLA allele as a “gene”. The resulting
expression matrix shows, for each cell, how many reads map to each allele — cells from
the PASC donor will express that donor’s private HLA alleles, and likewise for the
Healthy donor.
2.3.1 Preparing the HLA FASTA and GTF
The custom reference requires two inputs: a FASTA file of HLA allele sequences and a
matching GTF annotation file. The FASTA (HLA_Sequences_LK160824.fasta) contains one
entry per donor-specific HLA allele — each sequence is treated as a separate “chromosome”
in the custom reference, with the allele name as the sequence identifier.
The GTF is generated directly from the FASTA index rather than obtained from a database,
because the sequences are custom and not present in any standard annotation. samtools faidx first indexes the FASTA to produce a .fai file whose columns are sequence name,
length, and offset information. The awk command then reads the first two columns (name
and length) and writes a minimal GTF record for each allele, declaring a single exon
spanning the full length of the sequence. The gene_id and transcript_id attributes
are set to the allele name so Cell Ranger can link the GTF entries back to the FASTA
sequences.
Critical: Cell Ranger requires that no two features in the GTF have overlapping
genomic coordinates. Because all alleles are treated as features on their own separate
“chromosomes” this is not an issue as long as each allele has a unique sequence
identifier in the FASTA. If coordinate conflicts arise (e.g. two alleles sharing the
same sequence name), shift start positions incrementally (e.g. 1, 10000, 20000, …)
so that no two records overlap. Note that the second value in the .fai file is the
sequence length, not the end coordinate — the awk command uses it directly as the
GTF end field, which is correct.
# index the FASTA to generate the .fai file
# .fai columns: name, length, byte-offset, bases-per-line, bytes-per-line
samtools faidx HLA_Sequences_LK160824.fasta
# generate a minimal GTF from the FASTA index
# each allele gets one exon record spanning its full length (1 to length)
# gene_id and transcript_id are both set to the allele name
awk 'BEGIN {FS="\t"}; {
print $1 FS "gene" FS "exon" FS "1" FS $2 FS "." FS "+" FS "." FS \
"gene_id " $1 ";" "transcript_id " $1 ";"
}' HLA_Sequences_LK160824.fasta.fai > HLA_Sequences_LK160824.gtf2.3.2 Building the Cell Ranger reference and counting
# build a custom Cell Ranger reference from donor HLA sequences
~/Programs/cellranger-8.0.1/cellranger mkref \
--genome=HLA_sequences \
--fasta=HLA_Sequences_LK160824.fasta \
--genes=HLA_Sequences_LK160824.gtf
# quantify HLA allele expression per cell for each sample
~/Programs/cellranger-8.0.1/cellranger count --id=cellranger_sample_B_HLA \
--transcriptome=HLA_sequences \
--sample=cDNA-Sample-B-1 \
--fastqs=rawdata/ \
--localcores=16 \
--create-bam=true \
--chemistry=ARC-v1
~/Programs/cellranger-8.0.1/cellranger count --id=cellranger_sample_C_HLA \
--transcriptome=HLA_sequences \
--sample=cDNA-Sample-C-1 \
--fastqs=rawdata/ \
--localcores=16 \
--create-bam=true \
--chemistry=ARC-v1
~/Programs/cellranger-8.0.1/cellranger count --id=cellranger_sample_D_HLA \
--transcriptome=HLA_sequences \
--sample=cDNA-Sample-D-1 \
--fastqs=rawdata/ \
--localcores=16 \
--create-bam=true \
--chemistry=ARC-v12.4 HLA allele counting with scHLAcount
scHLAcount provides a second, complementary
approach to HLA-based donor confirmation. Rather than mapping all reads to a custom
reference, it directly counts molecules mapping to specific HLA alleles at the single-cell
level using a personalised reference genome derived from known HLA genotypes. It covers
class I genes (HLA-A, B, C) and class II genes (DPA1, DPB1, DRA1, DRB1, DQA1, DQB1).
The allele counts per cell are added to the Seurat object metadata and visualised per
SNP-assigned donor, providing an orthogonal confirmation that the donor0/donor1
labels correspond to the expected HLA profiles.
scHLAcount expects chromosome names without the chr prefix (i.e. 1, 2, … rather
than chr1, chr2, …), but the cellranger BAM files use UCSC-style names with the
prefix. The BAM files and barcode lists are therefore copied to a working directory and
the chromosome names are stripped with samtools and sed before running the tool.
The modified BAM files are then re-indexed with samtools index so scHLAcount can
perform random access by genomic position.
# working directory: scSNPdemux/scHLAcount/modify_cellranger/
# copy barcodes and BAM files from cellranger outputs into the working directory
cp ../../../cellranger_SC_ATAC_covid_B/outs/filtered_feature_bc_matrix/barcodes.tsv.gz cellranger_SC_ATAC_covid_B_barcodes.tsv.gz
cp ../../../cellranger_SC_ATAC_covid_C/outs/filtered_feature_bc_matrix/barcodes.tsv.gz cellranger_SC_ATAC_covid_C_barcodes.tsv.gz
cp ../../../cellranger_SC_ATAC_covid_D/outs/filtered_feature_bc_matrix/barcodes.tsv.gz cellranger_SC_ATAC_covid_D_barcodes.tsv.gz
cp ../../../cellranger_SC_ATAC_covid_B/outs/gex_possorted_bam.bam cellranger_SC_ATAC_covid_B_gex_possorted_bam.bam
cp ../../../cellranger_SC_ATAC_covid_C/outs/gex_possorted_bam.bam cellranger_SC_ATAC_covid_C_gex_possorted_bam.bam
cp ../../../cellranger_SC_ATAC_covid_D/outs/gex_possorted_bam.bam cellranger_SC_ATAC_covid_D_gex_possorted_bam.bam
# strip 'chr' prefix from all chromosome names, convert back to BAM, and re-index
# samtools view -h outputs SAM with header; sed removes the 'chr' prefix;
# the second samtools view -b converts back to binary BAM format
for i in *bam; do
samtools view -h $i | sed -e 's/chr//g' > ${i}_tmp.bam
samtools view -b ${i}_tmp.bam -o ${i}_mod.bam
samtools index ${i}_mod.bam
done
# run scHLAcount for each sample
# -d path to the HLA genotype database directory
# -o output directory for the allele count matrix
# --pl-tmp temporary directory for pseudoalignment intermediate files
./sc_hla_count_linux \
--bam modify_cellranger/cellranger_SC_ATAC_covid_B_gex_possorted_bam.bam_mod.bam \
--cell-barcodes modify_cellranger/cellranger_SC_ATAC_covid_B_barcodes.tsv.gz \
-d db_HLAgenotype/ -o hla_B --pl-tmp pseu_B
./sc_hla_count_linux \
--bam modify_cellranger/cellranger_SC_ATAC_covid_C_gex_possorted_bam.bam_mod.bam \
--cell-barcodes modify_cellranger/cellranger_SC_ATAC_covid_C_barcodes.tsv.gz \
-d db_HLAgenotype/ -o hla_C --pl-tmp pseu_C
./sc_hla_count_linux \
--bam modify_cellranger/cellranger_SC_ATAC_covid_D_gex_possorted_bam.bam_mod.bam \
--cell-barcodes modify_cellranger/cellranger_SC_ATAC_covid_D_barcodes.tsv.gz \
-d db_HLAgenotype/ -o hla_D --pl-tmp pseu_D