scgeno 0.1.5

Single-cell genotyping using Hidden Markov Model

scGeno - A Hidden Markov Model to denoise mouse genotype using Single-Cell RNA-seq data

What does scGeno do?

scGeno is a python package that uses a Hidden Markov Model to denoise the genotype segments of mammalian chromosomes (in this example, in mouse embryos), from single-cell RNA sequencing data of offspring organisms generated by crossing mixed genotype parental generation.

scGeno input

scGeno requires preprocessing of the single-cell data:

• Use SNPsplit to obtain unambiguously assignable reads for both genomes

SAMPLE=${1}         #/project/scRNA/processing/Exp
BC=${2}             #../outs/raw_gene_bc_matrices/mm10/barcodes.tsv(.gz)
R1=${3}             #"..._R1_001.fastq.gz ..._R1_001.fastq.gz ..._R1_001.fastq.gz"

R2=`echo ${R1} | sed 's/R1/R2/g'`

# get BCs surviving 10X pipeline
if file --mime-type ${BC} | grep -q gzip$; then
  zcat ${BC} | cut -f1 -d'-' >${SAMPLE}_10X_barcodes.tsv
else
  cut -f1 -d'-' ${BC} >${SAMPLE}_10X_barcodes.tsv
fi

# assignment of BC 2 read, reduce to 10X BCs
zcat ${R1} | sed 's/ .*//' | paste - - - - | cut -f1,2 | sed 's/^@//'  | perl -ane '$F[1]=substr($F[1], 0, 16); print "$F[0]\t$F[1]\n"' >${SAMPLE}_read.BC.tmp

fgrep -f ${SAMPLE}_10X_barcodes.tsv ${SAMPLE}_read.BC.tmp | sort >${SAMPLE}_read.BC.tsv
#rm ${SAMPLE}_read.BC.tmp

# join fastq files
cat ${R2} >${SAMPLE}_R2.fastq.gz

# alignment using STAR (parameters suggested by SNPsplit manual)
STAR --runThreadN 20 --genomeDir SNPsplit/SNPsplit_v0.3.2 --readFilesIn ${SAMPLE}_R2.fastq.gz --readFilesCommand zcat --alignEndsType EndToEnd --outSAMattributes NH HI NM MD --outSAMtype BAM Unsorted --outFileNamePrefix ${SAMPLE}_

# assignment of reads to genome1 (e.g. B6) or genome2 (e.g. CAST) using SNPsplit
perl SNPsplit --snp_file SNPsplit/SNPsplit_v0.3.2/all_SNPs_CAST_EiJ_GRCm38.txt.gz ${SAMPLE}_Aligned.out.bam --samtools_path samtools/samtools-1.6/

• Extract read IDs for genome1 (e.g. B6) or genome2 (e.g. CAST)

grep -w G1 sample1_unambiguous.bed | cut -f4 >sample1_G1.txt
grep -w G2 sample1_unambiguous.bed | cut -f4 >sample1_G2.txt

• Make folder

mkdir sample1_G1
mkdir sample1_G2

• Filter fastq I1/R1/R2/ for G1/G2

## sample_G1
seqkit grep --pattern-file sample1_G1.txt sample1_I1_001.fastq.gz >sample1_G1/sample1_G1_I1_001.fastq
seqkit grep --pattern-file sample1_G1.txt sample1_R1_001.fastq.gz >sample1_G1/sample1_G1_R1_001.fastq
seqkit grep --pattern-file sample1_G1.txt sample1_R2_001.fastq.gz >sample1_G1/sample1_G1_R2_001.fastq
gzip sample1_G1/*.fastq

## sample1_G2
seqkit grep --pattern-file sample1_G2.txt sample1_I1_001.fastq.gz >WT65_G2/WT65_G2_I1_001.fastq
seqkit grep --pattern-file sample1_G2.txt sample1_R1_001.fastq.gz >WT65_G2/WT65_G2_R1_001.fastq
seqkit grep --pattern-file sample1_G2.txt sample1_R2_001.fastq.gz >WT65_G2/WT65_G2_R2_001.fastq
gzip sample1_G2/*.fastq

• Run cell ranger count separately on G1 and G2

• Run cell ranger aggregate separately for G1 and G2 if needed

We can now install scGeno via pip:

pip install scGeno

After installation we can look at the help page and learn about input and output files:

scGeno --help

Now we can run scGeno! These are the inputs the program needs:

• aggregated feature barcode matrix, already filtered after standard scRNA-seq QC
• G1 feature barcode matrix
• G2 feature barcode matrix
• simplified gtf (csv) file of the specie we are interestes in, e.g.:

chr1,gene,3073253,3074322,"TEC","4933401J01Rik"
chr1,gene,3102016,3102125,"snRNA","Gm26206"
chr1,gene,3205901,3671498,"protein_coding","Xkr4"
chr1,gene,3252757,3253236,"processed_pseudogene","Gm18956"
...

This file for mm10 is located in the metadata folder.

• known imprinted genes list of the specie we are interested in An example of this file for mm10 is located in the metadata folder.
• the numbers of replicates and chromosomes you want to use to train the model

Example run:

scGeno -s WT_mm10 -gtf mm10.simple.gtf -imp mm10_imprinted_genes.csv -rep embryo -a WT_aggr.h5 -G1 WT_G1.h5ad -G2 WT_G2.h5ad -r 8 -c 9 -o out_path 

scGeno will output the following files:

• one dataframe containing for all the replicates per sample (columns), and for all genes (rows), 0 or 1 which indicates the genotype state. In the example reported here: 0 = mixed genotype, 1 = not mixed genotype (*_allreps.csv)
• one dataframe per replicate having these columns:

chr	gene	ratio	observed_state	hidden_state

• one bed file per replicate indicating the chromosomal segments having a continuous genotype state