comprehension for 'Using prephased SNP profiles'

[hanjun98]

Hi, first, Thank you for awesome tool to create new insight! I felt person who made numbat is so meticulous.

I have a question to 'Using prephased SNP profiles'

''' Using DNA-derived genotype information is another way to improve SNP density and phasing. If you have SNP calls from DNA genotyping (e.g. WGS/WES), you can first perform phasing on the DNA-derived VCF. Then run cellsnp-lite on scRNA-seq BAMs against the DNA-derived VCF to generate allele counts (only include heterozygous SNPs). Finally, merge the phased GT fields (from phased DNA-derived VCF) with the obtained allele counts to produce an allele dataframe in the format of df_allele (see section Preparing data). '''

1. If you have SNP calls from DNA genotyping (e.g. WGS/WES), you can first perform phasing on the DNA-derived VCF. Then run cellsnp-lite on scRNA-seq BAMs against the DNA-derived VCF to generate allele counts (only include heterozygous SNPs). --> I already have VCF(derived from WES) so I performed phasing on the DNA-derived VCF like this

`pileup_R='/mnt/gmi-l1/_90.User_Data/dhthxkr/codes/MDS/sc-seq/numbat/pileup_and_phase.R' gmap='/mnt/gmi-l1/_90.User_Data/dhthxkr/tools/Eagle_v2.4.1/tables/genetic_map_hg38_withX.txt.gz'

snpvcf='/mnt/gmi-l1/_90.User_Data/dhthxkr/ref/numbat/genome1K.phase3.SNP_AF5e2.chr1toX.hg38.vcf'

paneldir='/mnt/gmi-l1/_90.User_Data/dhthxkr/ref/numbat/1000G_hg38' eagle='/mnt/gmi-l1/_90.User_Data/dhthxkr/tools/Eagle_v2.4.1/eagle'

sample=$1 WES=$2

mkdir -p ${sample}/Numbat_WES

Rscript ${pileup_R} \ --label ${sample} \ --samples ${sample} \ --bams ${sample}_cellranger/outs/possorted_genome_bam.bam \ #single cell cellranger ouput bam --barcodes ${sample}/leiden_v2.4/matrix_files/barcodes.tsv.gz \ #single cell barcodes --outdir ${sample}/Numbat_WES \ --gmap ${gmap} \ --eagle ${eagle} \ --snpvcf ${sample}/WES/${WES}.vcf.gz \ #WES derived VCF(contains germlines and somtics, only filter='PASS') --paneldir ${paneldir} \ --ncores 5 &`

so I can get a '${sample}/Numbat_WES/MDSC02_allele_counts.tsv.gz'

2. Finally, merge the phased GT fields (from phased DNA-derived VCF) with the obtained allele counts to produce an allele dataframe in the format of df_allele. --> what's mean 'merge the phased GT fields'? Actually, I inferred that merging allele_counts.tsv(Pileup output with option(--snpvcf ${sample}/WES/${WES}.vcf.gz)) with allele_counts.tsv(Pileup output with option(--snpvcf genome1K.phase3.SNP_AF5e2.chr1toX.hg38.vcf))

so, I make a merged_allele_counts.tsv(from option(--snpvcf ${sample}/WES/${WES}.vcf.gz) allele_counts.tsv and option(--snpvcf genome1K.phase3.SNP_AF5e2.chr1toX.hg38.vcf)) allele_counts.tsv

And then run 'run_numbat' with merged_allele_counts.tsv like this

`#!/usr/bin/env Rscript args = commandArgs(trailingOnly=TRUE) sample <- args[1]

library(glue) library(stringr) library(data.table) library(dplyr) library(vcfR) library(Matrix) library(numbat) library(pagoda2) library(ggplot2) library(ggtree) library(tidygraph) library(patchwork)

df_allele=fread(paste0('/mnt/gmi-l1/_90.User_Data/dhthxkr/MDS/',sample,'/Numbat_WES_v2/',sample,'_merged_allele_counts.tsv.gz'), header = T) count_mat = readMM(paste0('/mnt/gmi-l1/_90.User_Data/dhthxkr/MDS/',sample,'/QC_output/matrix_files_v2.4/matrix.mtx.gz')) cells = fread(paste0('/mnt/gmi-l1/_90.User_Data/dhthxkr/MDS/',sample,'/QC_output/matrix_files_v2.4/barcodes.tsv.gz'), header = F)$V1 genes = fread(paste0('/mnt/gmi-l1/_90.User_Data/dhthxkr/MDS/',sample,'/QC_output/matrix_files_v2.4/features.tsv.gz'), header = F)$V2 colnames(count_mat) = cells rownames(count_mat) = genes count_mat = as.matrix(count_mat) count_mat = rowsum(count_mat, rownames(count_mat)) count_mat = as(count_mat, "dgCMatrix")

out = run_numbat( count_mat, # gene x cell integer UMI count matrix ref_hca, # reference expression profile, a gene x cell type normalized expression level matrix df_allele, # allele dataframe generated by pileup_and_phase script genome = "hg38", t = 1e-5, ncores = 20, plot = TRUE, out_dir = paste0('/mnt/gmi-l1/_90.User_Data/dhthxkr/MDS/',sample,'/Numbat_WES_v2') )`

But I got this error "Warning message: In asMethod(object) : sparse->dense coercion: allocating vector of size 2.4 GiB Error in check_allele_df(df_allele) : Inconsistent SNP genotypes; Are cells from two different individuals mixed together? Calls: run_numbat -> check_allele_df Execution halted"

So, I don't know when it wrong?

Thank you for reading! I'll be waiting for your answer!