Single-cell level haplotype phasing is key to studying clonal hematopoiesis, X chromosome inactivation, RNA editing, mitotic recombination, and somatic mutations. The availability and exponential growth of scRNA-seq data represent a vastly untapped resource for read-backed phasing in single cells. Unlike bulk RNA-seq, scRNA-seq’s read structure, high PCR duplication, and low sequencing depth present challenges. We introduce euler, a method that accurately infers long-range haplotypes by leveraging the read structure and PCR duplicates unique to scRNA-seq data. euler can be utilized in these following three scenarios:
- Germline variants, where all cells share the same haplotypes;
- Somatic variants, where haplotypes are shared by mutant cells and their progenies;
- RNA editing variants, where haplotypes are unique to each cell. euler uses ultrafast graph algorithms and is scalable to atlas-scale data with 1M single cells.
Developed and maintained by Laurentius.
Runs on python 3.x, requires samtools, vcftools, bedtools.
- euler - cEll-resolved Umi-Linked haplotypE Reconstruction
Overall Structure
The requirement is the same as installing from pypi.
$ git clone https://github.com/boxiangliulab/euler.git
$ cd euler
$ pip install .Note only install from source code is available before formal publication.
The typical installation time is ~1 min depends on your on computer/HPC.
# On all Chromosomes
$ euler --id <ID> --bam <path-to-bam> --vcf <path-to-vcf> --sample <SAMPLE_NAME> --total_chr <number-of-chrs-without-chrX> --input_type <cellranger/umitools/star/re>
# On one Chromosome
$ euler --id <ID> --bam <path-to-bam> --vcf <path-to-vcf> --sample <SAMPLE_NAME> --chr <desired_chr> --input_type <cellranger/umitools/star/re> Requires a VCF and BAM, produces a HAPCUT-style output and a phased VCF with computed haplotype phases and result files containing haplotype details, graph connectivities, and read counts.
We highly recommend you to check your environment setup using the provided example data under example/ directory.
Run command
$ cd /path/to/euler
$ euler --id test --bam ./example/bam/NA12878.chr5.example.bam --vcf ./example/vcf/NA12878.example.vcf.gz --sample NA12878 --chr chr5 --input_type umitools If everything was set up right, this command should run with out any error.
For more example inputs and example inputs for euler-somatic, refer to here. After downloaded the example input files, try out the given bash scripts to validate the installation of euler.
euler supports BAM files generated by various single-cell alignment/analysing tools e.g. cellranger, star-solo, umitools etc, and multiple filtering parameters on BAM file including alignment score and MAPQ.
Apart from inputing single BAM file using --bam <path-to-bam>, we provided a way to input multiple BAM files at once. This will be helpful when mixing different reads from tissures of a common individual, e.g. normal cells and tumor cells.
First, you need to prepare a bam.list file
#bam.list
<path-to-bam1>,ID1
<path-to-bam2>,ID2
...and use the argument --bam_list bam.list. The ID column will be append to cell barcodes from corresponding BAM file.
We defined cell barcode and UMI style of different softwares as follows:
#Cellranger
SRR8551677.314853539 ... CB:Z:<Barcodes> UB:Z:<UMI> ...
#UMI-tools
SRR8551677.236615618_<Barcodes>_<UMI> ...
#star-solo
<Barcodes>_<UMI>_AAAAAEEEEE_SRR6750053.36514992 ...You shall specify your input type by using --input-type <cellranger/umitools/star>.
Noted that the cell barcode and UMI style may change between different versions of the same software, and there may be more styles out of our selected ones, we present you with --input-type re, which allows you to use custom regular expressions to extract UMI and cell barcodes.
For example, given a BAM file, one line of the BAM file show as
SRR8551677.64678072 16 chr1 14426 1 ... BC:ACTACTACT UMI:TAGTAGTGA ...
The corresponding regular expression for extract cell barcode and UMI from this BAM is
--bc_re "BC:([AGCT]+)" --umi_re "UMI:([AGCT]+)"
Further information of regular expression can be found here.
We applied multiple threshold on the BAM file to achieve balance between precision and sensitivity, with corresponding argument --mapq, --as_quality etc. Detailed information can be seen here
Variants pending for phasing are given by VCF files. By providing VCF files using --vcf, you should provide the desired sample name in the VCF file with --sample.
Sometimes there are no ground truth VCF files, nor is there an existing VCF file for an individual. In other words, you only have the scRNA-seq bam file. Worry not! You can call variants on your own using the following command.
$ samtools index <path-to-bam>
$ bcftools mpileup -Ou -f <path-to-genome.fa> <path-to-bam> | bcftools call -mv -Ob -o tmp.bcf
$ bcftools view tmp.bcf -o tmp.vcf
$ bgzip -c tmp.vcf > <your-desired-name>.vcf.gz
$ tabix -p vcf <your-desired-name>.vcf.gzYou can perform further filtering on this VCF file or simply indicate --raw_vcf and --vcf_qual <threshold-on-vcf-QUAL> when using euler.
Note that samtools and bcftools are not optimized for single-cell data. Consider using Monopogen, which designed for scRNA-seq data. A typical output of Monopogen is like:
├── germline
│ ├── chr1.gl.vcf.gz
│ ├── chr1.gp.vcf.gz
│ ├── chr1.phased.vcf.gz
├── Script
└── somatic
├── chr1.allSNVs.csv
├── chr1.cell_snv.cellID.csv
├── chr1.cell_snv.cellID.filter.csv
├── chr1.cell_snv.mat.gz
├── chr1.cell_snv.snvID.csv
├── chr1.germlineTrioLoci_model.csv
├── chr1.germlineTwoLoci_model.csv
├── chr1.gl.filter.hc.cell.mat.gz
├── chr1.gl.vcf.DP4
├── chr1.gl.vcf.filter.DP4
├── chr1.gl.vcf.filter.hc.bed
├── chr1.gl.vcf.filter.hc.pos
├── chr1.putativeSNVs.csv
├── chr1.SNV_mat.RDS
├── LDrefinement_germline.chr1.pdf
└── svm_feature.chr1.pdf
A typical Monopogen output file chr<X>.putativeSNVs.csv may look like:
chr,pos,Ref_allele,Alt_allele,Depth_total,Depth_ref,Depth_alt,SVM_pos_score,LDrefine_twoLoci_score,LDrefine_trioLoci_score,LDrefine_merged_score,BAF_alt
chr1,633887,G,A,329,2,327,0.660488303686431,NA,NA,NA,0.990909090909091
chr1,634012,T,C,2402,2388,14,0.173676949268894,NA,NA,NA,0.00582605076987099
chr1,634133,C,T,431,401,30,0.25421399379009,NA,NA,NA,0.0694444444444444
chr1,958944,C,T,80,79,1,0.235030370462956,NA,NA,NA,0.0123456790123457To use Monopogen's output as euler's input, one should indicate --mono. An example command is
$ tabix -p vcf <path-to-monopogen-output>/germline/chr1.gp.vcf.gz
$ euler --id <ID> --bam <path-to-bam> --vcf <path-to-monopogen-output>/germline/<chr>.gp.vcf.gz --mono <path-to-monopogen-output>/somatic/<chr>.putativeSNVs.csv --sample <SAMPLE_NAME> --chr <desired_chr> --input_type <cellranger/umitools/star/re> If --mono is specified, somatic variants and germline variants will be annotated both in .euler output and output .vcf files.
Somatic variants called by Monopogen is filtered in euler. euler performs filtering on Depth_ref with --mono_ref, Depth_cal with --mono_alt, SVM_pos_score with --mono_svm, LDrefine_merged_score with --mono_ld and BAF_alt with --mono_baf_u and --mono_baf_l.
Sometimes, the input chromosome may contain the "_" character, e.g. mm10_1. To prevent errors, you need to specify --sep <any-character-you-like-except-from-_-and-:>. And in other circumstances, the naming style of given BAM file and VCF file may be different. You can use --chr_vcf <chr-on-VCF-file> to specify the desired chromosome in the VCF file.
euler assumes diploid for input, thus if there are more than one ALT for a input variant, the default parameters would lead to misleading results. If you consider this as unacceptable, or difference between ALT alleles are important in your study, set --non_binary_variant. Therefore the euler is provided ability to distinguish between different ALT.
Typical output consists of a single chr*.euler file. The format of the chr*.euler is in consistent with HAPCUT and HapTreeX. A typical block of the chr*.euler is presented as follow
...
BLOCK: offset: 249337 len: 2 phased: 2 SPAN: 75 correct: 1
0 1 0 chr20 249337 C G 0|1
1 0 1 chr20 249412 C T 1|0
****************From left to right, each corresponding column represents umique_id, allele on one haplotype, allele on the other haplotype, chromosome, position, REF, ALT and GT from the input VCF file.
If you would like to output the phased vcf file, you can specify --output_vcf. For the output vcf, we added the following fields to phased variants.
- SG: euler Local Genotype
- SB: euler Local Block
- PS: Phased Set Index (unique for each block)
A typical output VCF file would be like:
...
chr1 1023775 . G A 0 PASS KM=9.62;KFP=0;KFF=0;MTD=bwa_freebayes,bwa_gatk,bwa_platypus,isaac_strelka GT:SG:SB:SI 0|1:1|0:chr1_1023775_._G_A;chr1_1023789_._G_C;chr1_1023851_._G_A:47820887733280
chr1 1023789 . G C 0 PASS KM=9.61;KFP=0;KFF=0;MTD=bwa_freebayes,bwa_gatk,bwa_platypus,isaac_strelka GT:SG:SB:SI 0|1:1|0:chr1_1023775_._G_A;chr1_1023789_._G_C;chr1_1023851_._G_A:47820887733280
...
If the output of the ALG and the edge information of it is required, specify --allele_linkage. The output of edge linkage information between variants of each cell is outputed in <output_dir>/<id>/cell_allele_connections_chr*.txt. The typical content of the output file is
barcode var geno support
CCACGGAAGACAATAC chr22_16614653_._C_T,chr22_16614872_._A_T 00 1
TGAGCCGGTTCCCGAG chr22_16614653_._C_T,chr22_16614872_._A_T 00 1
ACTTACTGTTCCTCCA chr22_16614653_._C_T,chr22_16614872_._A_T 00 1
...
To enable residual graph output, specify --residual_edges. The results are stored in <output_dir>/<id>/singular_cells/singular_cell_linkage_*.txt. The typical content would be like
barcode var geno support p-value oppo_support total_covered_barcodes global_oppo_support correct is_singular
CCATTCGCAGGATTGG chr15_22690417_._A_G,chr15_22690362_._A_G 01 1 0.8607110296641832 0 7 5 1 0
CTCAGAAAGGAATTAC chr15_22690417_._A_G,chr15_22690362_._A_G 01 1 0.8607110296641832 0 7 5 1 0
where support represents the number of mRNAs supporting this pair of allele in this specified cell barcode. oppo_support represents the number of mRNAs in conflict with this pair in this cell barcode. total_covered_barcodes represents the number of cells that expressed this two variants. global_oppo_support represents the number of mRNAs in conflict with this pair in all cells.
To perform further somatic mutation analysis with euler-somatic, the user shall indicate --output_conflict when calling euler. The output conflict graph is stored in <output_dir>/<id>/conflicted_graphs/. Each .pickle file contains a networkx.Graph() object, user can easily load the graph using pickle.load() and perform custom analysis.
For detailed output, and graphical demonstration of ALGs (allele linkage graph), specify --verbose. All the conflicted ALGs will be stored in <output_dir>/<id>/conflict_graphs_graphml/*.graphml. The resolved conflicted ALGs will be stored in <output_dir>/<id>/resolved_conflict_graphs_graphml/*.graphml. The .graphml file can be further visualized and analysed using Gephi.
Although Monopogen has provided filtering strategy for filtering somatic mutation callset, there are still quite many sequencing errors hidden in the callset. In euler-somatic, we can utilize the phasing information to refine the somatic mutation callset, as sequencing errors may evenly distributed on both alleles.
To perform somatic mutation filtering, call euler-somatic filter and indicate the --id and --output_dir if the output directory is not ./output. A minimal example is
$ euler-somatic filter --id <id>This function will generate one output:
refined.somatic.mutations.tsv
An example output is shown below
chr19_45523951_somatic_T_A,somatic_mutation
chr12_7920572_somatic_A_G,sequencing_error
chr19_18589308_somatic_A_G,sequencing_error
chr19_8413314_somatic_T_G,sequencing_error
chr19_34860631_somatic_G_C,sequencing_error
chr12_7920565_somatic_A_C,sequencing_error
chr12_7934570_somatic_G_T,sequencing_error
chr19_40799190_somatic_C_G,sequencing_error
chr1_198658531_somatic_T_C,sequencing_error
chr1_1512167_somatic_C_G,sequencing_errorTo perform somatic variations analysis in the paper, you first need to specify --output_conflict when performing euler phasing.
We recommend users preprocess the GTF file and build index for genes with
$ euler-somatic init --gtf <path-to-gtf>To perform in-phase and out-of-phase detection, user shall first prepare a gtf file and run the following command.
# For somatic variation analysis
$ euler-somatic run --id <id> --gtf <path-to-gtf> --vcf <path-to-vcf> --chr <chromosome>An example GTF file could be downloaded here.
The script will generate two outputs:
chr*.in.phase.hits.annotated.csvchr*.out.of.phase.hits.annotated.csv
A typical example of output is
sample germline somatic gene cells supports
JP_RIK_H012 chr1_226070308_._C_T chr1_226067274_somatic_T_C H3-3A CGGTTAAAGCACCGTC-1_1 1
JP_RIK_H012 chr1_46676803_._A_G chr1_46677049_somatic_C_T EFCAB14 AGCTCTCTCATAGCAC-1_1 1
JP_RIK_H012 chr1_46706337_._C_T chr1_46677049_somatic_C_T EFCAB14 AGCTCTCTCATAGCAC-1_1 1
JP_RIK_H012 chr1_46718867_._C_T chr1_46677049_somatic_C_T EFCAB14 AGCTCTCTCATAGCAC-1_1 1
JP_RIK_H012 chr1_46718524_._A_G chr1_46677049_somatic_C_T EFCAB14 AGCTCTCTCATAGCAC-1_1 1
JP_RIK_H012 chr1_149402084_._C_CT chr1_149471694_somatic_A_G NOTCH2NLC GGGAATGTCTCTAAGG-1_1 2where each line represents a in/out-of-phase variant pair. germline column represents germline variants, somatic column represents somatic variants and gene column represents the gene name on which the event occurs. cells represents the cell barcodes that harbors the somatic-germline variant pair, and support represents the total cell counts.
Notablly, all the germline variants and somatic variants are not filtered by CADD score nor missense/synonymous. User may perform further filtering theirselves.
To perform somatic variations analysis in the paper, you first need to specify --allele_linkage.
If --allele_linkage is specified, euler will output variant pairs and its configuation in each cell barcode in <output_dir>/<id>/cell_allele_connections_<chr>.txt. A example is shown below
barcode var geno support
TGAAAGATCTGCTGCT chr22_17113996_._T_C,chr22_17114025_._G_A 00 1
GTCATTTGTGTCCTCT chr22_17113996_._T_C,chr22_17114025_._G_A 00 1
TCAGGATGTACATCCA chr22_17115288_._T_C,chr22_17115339_._T_C 11 1
CTCATTATCAACACCA chr22_17115288_._T_C,chr22_17115339_._T_C 11 1where barcode represents a unique cell, var represents a variant pair, geno represents the configuration of the given variant pairs (0 represents wild type and 1 represents alternative) and support represents the UMI count of the given variant configuration in the given barcode.
After allele linkage information is extracted, we can then perform RNA-editing analysis.
First, edit sites and potential edQTL sites shall be filtered. We perform this filtering by refering to the GTEx edQTL data. The GTEx edQTL data can be downloaded here.
After data preparation and filtered var into one RNA-editing site and potential edQTL, you shall count the summation of support with different geno and same barcode and var, which represents mRNA count of different configuration of edQTL and editing site pairs.
For example, assume in a certain cell_allele_connections_<chr>.txt
barcode var geno support
TGAAAGATCTGCTGCT <edQTL>,<Editing Site> 00 a
TGAAAGATCTGCTGCT <edQTL>,<Editing Site> 01 b
TGAAAGATCTGCTGCT <edQTL>,<Editing Site> 10 c
TGAAAGATCTGCTGCT <edQTL>,<Editing Site> 11 dthose four lines provide us the counting information of a certain pair of edQTL and editing site in a certain cell. We can calculate the editing level of <edQTL> on <Editing Site> by
ref_level = b / a + b
alt_level = d / c + d
Furthermore, we can generalize this into a certain cohort by group by the same var among all samples and analyze cell type-specific RNA-editing by annotating each barcode with its cell type and group by the annotated cell type. However, as data format and analysis aim may vary drastically for RNA-Editing, we could not provide a standard and generalized analysis pipeline for edQTL analysis.
The R script used to reproduce the RNA-Editing analysis in the euler paper is in example/euler_rna_editing.R for reference and further modifications.
We welcome any enquiries on usage of euler on single-cell edQTL analysis!
--id: A unique run ID string (e.g. NA12878).--bam: Indexed BAM file containing aligned reads.--bam_list: A list of input BAM files.--vcf: Input VCF file, must be gzipped and tabix indexed.--sample: Sample name in VCF.--chr: To restrict phasing in a given chr on BAM & VCF.--total_chr: Total chromosome count (without 'X' chromosome) for phasing on all chromosomes--raw_vcf: If the vcf is not filtered, and the FILTER column is not 'PASS', specify this option--sep: Character used to construct split variant information. Default_.--non_binary_variant: If set, means there may be multiple ALT. for a single variant in the VCF file. Defaultfalse.--chr_prefix: Chromosome prefix. Defaultchr.--chr_vcf: To restrict phasing in a given chr on VCF, if chromosome is not named equally between BAM and VCF--neglect_hla: If set, variants in the HLA region is neglected.--seed: Random seed.--tmp_dir: Directory of tempfile.--output_dir: Path to output directory. Output of euler will be stored in<output_dir>/<id>--input_type: How umi-barcode is provided, e.g. cellranger-style, umitools-style or 're' for custom regular expression.--bc_re: The regular expression for extracting Cell Barcode in the BAM file.--umi_re: The regular expression for extracting UMI in the BAM file.
--mapq_threshold: Threshold on minimum MAPQ value for input BAM. Reads with MAPQ lower than this value will be ignored.--vcf_qual: The quality threshold on QUAL during processing raw input VCF file.--interval_threshold: Edge with allele distance on genome more than this threshold will be ignored.--base_conflict_threshold: Determine whether different bp at the same position will be persisted.--method: Specify graph split method, e.g. mincut, fiedler.--fiedler_threshold: Nodes with corresponding value in fiedler vector lower than threshold will be removed--shortest_path: Decide whether activate split_graph_by_common_shortest_path..--remove_node: Activate when--shortest_pathis set. Remove no more thanremove_nodein split_graph_by_common_shortest_path()--as_quality: A filter on alignment score in BAM files.--edge_threshold: A filter on low confidence edges on graph.--layers: Number of GNN layers. Default 1. No GNN is applied if--layersis set to 0;
--verbose: Determine whether output conflicted graphs in pdf and graphml format.--benchmark: Determine wheter output benchmark information.--output_conflict: Decide whether to output conflict graphs. This is a must if you want to useeuler-somaticfor further somatic mutation analysis.--residual_edge: Decide whether to output residual edges (edges that are removed when resolving graph conflicts).--allele_linkage: Decide whether output allele linkage count.--output_vcf: Decide whether output vcf or not (severe performance decrease).
--thread: Number of thread number for multiprocessing. Default8.
--mono: The result putativeSNVs.csv of Monopogen somatic output, used for somatic variation analysis--mono_svm: Threshold on SVM_pos_score in putativeSNVs.csv. Default:0.5--mono_ld: Threshold on LDrefine_merged_score in putativeSNVs.csv. Default:0.5--mono_baf_u: The upper bound on BAF_alt in putativeSNVs.csv. Default:0.5--mono_baf_l: The lower bound on BAF_alt in putativeSNVs.csv. Default:0.1--mono_ref:Threshold on Dep_ref in putativeSNVs.csv. Default:5--mono_alt: Threshold on Dep_alt in putativeSNVs.csv. Default:5
--id: The ID ofeulerrun.--vcf: The germline VCF file of sampleid.--chr: Indicate the chromosome on which theeuler-somaticwill be performed on.--output_dir: If other--output_diris set when callingeuler, this should be set as the same--output_diraseuler.--gtf: Input GTF file. We recommend you to input a .gtf file rather than a .gtf.gz file, because euler-somatic will have to depress the file everytime you input a compressed gtf file, which will result in a waste of time and disk space.--cells: Threshold on number of supporting cells.--reads: Threshold on number of supporting reads.
--gtf: Input GTF file. We recommend you to input a .gtf file rather than a .gtf.gz file, because euler-somatic will have to depress the file everytime you input a compressed gtf file, which will result in a waste of time and disk space.
--id: The ID ofeulerrun.--output_dir: If other--output_diris set when callingeuler, this should be set as the same--output_diraseuler.
Not published yet. 😭