scgeno 0.1.5
Single-cell genotyping using Hidden Markov Model
License: MIT
Requires Python: >=3.10

# scGeno - A Hidden Markov Model to denoise mouse genotype using Single-Cell RNA-seq data
<div align="center">
  <img src="./images/Figure1_v1.svg" alt="High quality figure" width="70%">
</div>

### 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