Watch Out for a Second SNP: Focus on Multi-Nucleotide Variants in Coding Regions and Rescued Stop-Gained

Abstract

Most single-nucleotide polymorphisms (SNPs) are located in non-coding regions, but the fraction usually studied is harbored in protein-coding regions because potential impacts on proteins are relatively easy to predict by popular tools such as the Variant Effect Predictor. These tools annotate variants independently without considering the potential effect of grouped or haplotypic variations, often called “multi-nucleotide variants” (MNVs). Here, we used a large RNA-seq dataset to survey MNVs, comprising 382 chicken samples originating from 11 populations analyzed in the companion paper in which 9.5M SNPs— including 3.3M SNPs with reliable genotypes—were detected. We focused our study on in-codon MNVs and evaluate their potential mis-annotation. Using GATK HaplotypeCaller read-based phasing results, we identified 2,965 MNVs observed in at least five individuals located in 1,792 genes. We found 41.1% of them showing a novel impact when compared to the effect of their constituent SNPs analyzed separately. The biggest impact variation flux concerns the originally annotated stop-gained consequences, for which around 95% were rescued; this flux is followed by the missense consequences for which 37% were reannotated with a different amino acid. We then present in more depth the rescued stop-gained MNVs and give an illustration in the SLC27A4 gene. As previously shown in human datasets, our results in chicken demonstrate the value of haplotype-aware variant annotation, and the interest to consider MNVs in the coding region, particularly when searching for severe functional consequence such as stop-gained variants.

Introduction

Next-generation sequencing has given access to genomes at the nucleotide level through DNA-seq but also specifically to expressed regions by whole-exome sequencing (WES, originally focusing on exonic parts of the genome) or RNA-seq. These data enable us to call genetic variations by spotting differences between aligned reads and the species reference genome or among aligned reads. Among these genetic variations, single-nucleotide polymorphisms (SNPs) are the most frequent and most studied variations. Although most variations are located in non-coding regions, the most analyzed lie in protein-coding regions where their potential impact(s) on the protein are relatively easy to predict. For example, in a study using 60,706 human exomes, the Exome Aggregation Consortium (ExAC) identified 3,230 genes with near-complete depletion of predicted protein-truncating variants. Of these genes, 72% have not been related to any known human disease phenotype (Lek et al., 2016). Different popular tools have been developed this last decade to predict SNPs’ effects on proteins such as Variant Effect Predictor (VEP) (McLaren et al., 2016), SnpEff (Cingolani et al., 2012), or ANNOtate VARiation (ANNOVAR) (Wang et al., 2010). But these tools consider each variation location individually, as if it they were specific to “reference” nucleotides. However, SNPs can be grouped by two or more coexisting variants present in the same haplotype (in the same individual), in which case they are called “multi-nucleotide variants” (MNVs). An example of MNV in one individual (with two nearby SNPs) is given in Figure 1. When such MNVs occur within a codon, the amino acid modification caused by this MNV may be different from protein change resulting from each constituent SNPs taken individually, leading to a risk of erroneous functional consequence prediction, as depicted in Figure 1. MNV identification tools have been developed using different methods for phasing SNPs [MAC (Wei et al., 2015), varDic (Lai et al., 2016), COPE (Cheng et al., 2017), BCFtools (Danecek and McCarthy, 2017), and MACARON (Khan et al., 2018)] and have been applied to different human genetic variant datasets [1,000 Genomes Project dataset (Cheng et al., 2017; Danecek and McCarthy, 2017; Khan et al., 2018; Wang et al., 2020), ExAC (Lek et al., 2016), The Cancer Genome Atlas (Lai et al., 2016), or gnomAD consortium (Wang et al., 2020)], mainly based on exomes.

FIGURE 1

Robinson et al., 2011) screen shot indicates the principle of read-based phasing of SNPs: short read mapping against the reference genome of the heterozygous individual allows us to phase both SNPs giving two haplotypes: C with T (reference alleles) on one side and T with A (alternative alleles) on the other side. When translated, these two haplotypes correspond to a leucine or a STOP codon and not to a simple amino acid change (LEU → GLU) if the two haplotypes had been composed by only one reference and one variant as shown in (A). (C) Information (PID and PGT) provided by GATK in the VCF files about the phased SNPs according to the read-based phasing shown in (B).

To our knowledge, no study has been conducted on MNVs in livestock species. The aim of this paper is to focus on MNVs occurring in protein-coding regions to provide examples and evaluate the functional consequences of resulting mis-annotations. Considering this aim, we used 9.5M SNPs recently detected in 382 chickens from 767 multi-tissue RNA-seq, enriched by construction in expressed regions and therefore in protein-coding regions. From this 9.5M SNPs, we focused on the 3.3M SNPs with reliable genotypes [see the companion paper (Jehl et al., 2021)]. MNV identification requires properly phased variants, i.e., to be located either on the same haplotype (called therefore MNV) or on two different haplotypes (a case of individual SNPs) (see Figures 1A,B). Different SNP phasing strategies exist: (i) population-based phasing, using statistical inference of phase from haplotypes shared among individuals of a large genotyped population; (ii) family-based phasing, which analyzes the co-transmission of variants between parents and offsprings; and (iii) read-based phasing, which evaluates whether close variants are present on the same reads in the DNA-seq or RNA-seq data. Read-based phasing is particularly relevant for close variants, making this method appropriate for MNV analysis in codons, in which variants fall within a maximum distance of 2 bp from one another. Therefore, in this study, we have chosen to identify MNVs by using the read-based phasing provided by the HaplotypeCaller tool of the GATK toolkit, in the VCF file through additional fields (Figure 1C, with the PID and PGT fields) recently added by the common variant caller (DePristo et al., 2011; Van der Auwera et al., 2013; McKenna et al., 2020).

Materials and Methods

SNP Dataset

The 3,276,615 SNPs analyzed in this study have been detected following the method presented in the companion paper (Jehl et al., 2021) using 767 multi-tissue RNA-seq of 382 birds from 11 chicken populations (see Additional File 1). This SNP set corresponds to the union of the SNPs with reliable genotypes found in each population (list available on http://www.fragencode.org/lnchickenatlas.html). Briefly, variant detection was performed for each sample using the HaplotypeCaller tool of GATK toolkit (DePristo et al., 2011; Van der Auwera et al., 2013; McKenna et al., 2020) 3.7.0 with options “--stand_call_conf 20.0,” “--min_base_quality_score 10,” and “--min_mapping_quality_score 20” (which are the defaults values). The “GenotypeGVCFs” function was then used with the option “--stand_call_conf 20.0,” to jointly genotype all these samples into one VCF per tissue. First, biallelic SNPs were then extracted using the “SelectVariant” function with option “--selectType SNP –restrictAllelesTo BIALLELIC.” Variants were filtered using “VariantFiltration” with “QD < 2” and “FS > 30.” Considering genotypes, variants were selected with a “(5.reads.DP) genotype CR ≥ 20%” and a “CR ≥ 50%.” The 11 populations include a red jungle fowl population (RJFh), three broiler populations with one experimental line (FLLL) and two commercial ones (Cobb, HerX), six layer populations with two brown-egg commercial lines (Novo1 and Novo2), two brown-egg experimental populations (RpRm and LSnu), two white-egg or cream-egg experimental populations (FrAg and FAyo), and finally a cross between white- and brown-egg experimental lines (Rmx6).

Analysis of the Functional Impact of Each Individual SNP in the Coding Regions

VEP v92 (McLaren et al., 2016) with a GTF file enriched in long non-coding genes (“--gtf”) was used for effect prediction of each SNP, with “--everything” and “--total_length” options to respectively, obtain SIFT score predictions and lengths of cDNA, CDS, and protein positions (Ng and Henikoff, 2003; Sim et al., 2012).

MNV Calling and Recalculation of Consequences

The script to detect the MNV and to calculate the consequences is available in Additional File 2.

Detection of SNPs Located in the Same Codon

With the information produced by VEP, an ID composed of the “transcriptID” and the “position of the SNP in the coding sequence (expressed in codon number)” was created for each consequence. Through this approach, the same codon of the same transcript supporting at least two different SNPs will have the same ID. Thus, only duplicated IDs were kept as they correspond to those containing two or more SNPs.

Detection of Co-located SNPs Carried by the Same Haplotype (MNV)

To test if the SNPs located in the same codon were also present in the same haplotype, the VEP file generated in the previous step and the VCF file were joined by the “SNPid” key, equivalent to “CHR_POS_REF/ALT.” The resulting file (VEP merged to VCF information) contained SNPs on the same codon with additional information about their phase (PID and PGT). Finally, the SNPs which were phased (i.e., same PID) and co-located in a codon were extracted: they correspond to MNVs containing two or three phased SNPs in the same codon.

Recalculation of the Consequences

With the R package Biostrings v2.50.2 (Pagès et al., 2021), the associated amino acids were produced for each MNV, and with the same strategy as VEP being adopted, the MNV consequences were established.

Analysis of MNV Functional Impacts and Comparison With the Constituent SNP Impacts

To compare MNV and independent SNP consequences, we selected only the most impactful consequence per codon for these constituent SNPs using the following order of priority from severe to weak consequences: (1) stop-gained, (2) stop-lost, (3) start-lost, (4) missense_variant, (5) stop-retained_variant, and (6) synonymous_variant.

For MNVs with a missense annotation corresponding to a missense annotation for both constituent SNPs, we distinguished two cases:

• –

  missense MNV with an amino acid different from those predicted by the constituent SNPs (SNP1: Missense A; SNP2: Missense B → MNV: Missense C) and

• –

  missense MNV with an amino acid common to one of two amino acids predicted by the constituent SNPs (SNP1: Missense A; SNP2: Missense B → MNV: Missense A or B).

In order to visualize the results, we produced an alluvial plot using the R “alluvial” package v0.1-2 (Bojanowski, 2020).

GO or KEGG Term Enrichment Analysis

The enrichment analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) terms in the gene set of interest was performed using the STRING v11.0. tool (Szklarczyk et al., 2019), and a GO or KEGG term was found significantly enriched if the BH-adjusted p ≤ 5%.

DNA Sequencing of SLC27A4

Five microliters of DNA samples was mixed with 5 μl of GoTaq Flexi Buffer 5 ×, 2 μl of MgCl₂ solution (25 mM), 0.125 μl of GoTaq DNA polymerase (5 U/μl) (Promega, catalog number: M891), 0.5 μl of dNTPs 10 mM, 12.5 μl H₂O, and 1.25 μl of specific reverse (CATTCCCGTAGTGCCAGAGG) and forward primers (GCACTTTCTGGTGCAAAGCA) at 10 μM. Reaction mixtures were then incubated in a T100 thermal cycler (Bio-Rad, Marne la Coquette, France) for 30 cycles with 30 s at 94°C, 30 s at 60°C, and 30 s at 72°C. The amplification products were then deposited on a 2% agarose gel and sent for sequencing (Genoscreen) to verify their location on the chicken genome.

Results and Discussion

Read-Based Phasing for Identification of MNVs

Using 3.3M SNPs previously detected from 767 multi-tissue RNA-seq of 382 animals from 11 chicken populations and therefore enriched in coding regions [see the companion paper (Jehl et al., 2021), section “Materials and Methods”], we identified 260,919 unique SNPs in 26,702 transcripts corresponding to 15,835 genes out of 19,545 protein-coding genes (Figure 2, right part—in yellow).

FIGURE 2

Jehl et al., 2021). Right: MNV constituent SNP selection and protein impact selection of these SNPs separately analyzed by VEP.

As shown in Figure 2 (left part—in green), we then defined an MNV in a codon as a group of two or three phased SNPs, i.e., existing on the same haplotype in the same individual. We found 11,183 SNPs (4.3% of the SNPs in codons) as constituent variants of 5,533 MNVs, which corresponded to 4,415 transcripts and 2,916 genes. Most of them (98%: 5,416) contained two SNPs with similar proportions (1/3) by constituent SNP position in the codon (1–3, 1–2, and 2–3, Figure 2, left). In order to ensure the reliability of the MNVs, we selected MNVs observed in at least five individuals. Out of the 5,416 MNVs with two SNPs, 2,965 MNVs were present in at least five individuals, corresponding to 2,636 transcripts and 1,792 genes. No GO terms or KEGG terms were found as significantly enriched for this gene list, suggesting that no specific biological pathway was impacted by MNVs. Table 1 gives the distribution of MNVs and their consequences according to the individual number supporting the MNV (ranging from 2 to 100 individuals). We can note that 31% of the 5,416 MNVs with two SNPs are observed only in a single individual and are here considered as erroneous, likely due to sequencing errors.

Functional Impact Comparison of MNVs and Their Component SNPs

Focusing on the 2,965 MNVs present in at least five individuals, we then compared their functional consequences with those of the 5,930 constituent SNPs as illustrated by the right part of the workflow provided in Figure 2. For such a comparison, we retained for each MNVs, the most severe consequence of the constituent SNPs according to the order indicated in Figure 2 (bottom right). The alluvial plot in Figure 3 depicts the consequence variations before (left) and after (right) taking the MNV impacts into account, according to the different consequence categories; the details of impact variation per MNV are given in Additional File 3 for the whole 2,965 MNV set. We can observe in Figure 3 that the biggest change in variant impacts concerns the originally stop-gained consequence categories, for which 95.6% were re-predicted as missense (green flux: 87 out of the 91 stop-gained initially predicted). The second and third biggest fluxes concern missense consequence categories, for which 37.3% had a different predicted amino acid (violet flux: 1,038 MNVs out of the 2,780 initial missenses), and 3.0% became synonymous variants (blue flux: 83 MNVs out of the initial missenses). The distribution of re-prediction fluxes is provided in Table 1 as a function of the individual number supporting the MNV among the 382 individuals analyzed. Among the 87 rescued stop-gained observed in five individuals, half (47) are observed in at least 15 individuals and are present on average in five populations (see Additional File 4). Out of the MNVs, the proportion of rescued stop-gained MNVs (2.9%), defined as at least one of the individual SNPs creating a nonsense mutation but not the resulting MNV, is in the same order of magnitude as the one reported by the gnomAD consortium with 1,821 rescued stop-gained MNVs out of 31,575 human MNVs (5.8%) (Wang et al., 2020). Genes with a stop-gained MNV rescued in the missense variant are available in Additional File 4 with the population affected and the individual number per population carrying these MNVs. To a lesser extent, nine missenses were re-predicted as stop-gained, which would have gone unnoticed without re-prediction. After a deeper investigation with the IGV browser, these re-predicted stop-gained variants seem to be present since they were not located in a potential exon skipping. Finally, this stop-gained category drastically declined by 86% (from 91 to 13) after considering MNVs, whereas the synonymous category was increased by twofold (from 79 to 159). These different category changes after considering MNVs have a major impact on variant interpretation and thus are critical for accurate variant annotation. More broadly, when the MNVs were considered together, the resulting functional impact differed from the independent impacts of the individual variants in 41.1% of the analyzed MNVs. This large percentage of mis-annotations is relatively consistent with ∼60% of reannotations in human MNVs recently reported by the gnomAD consortium in coding regions (Wang et al., 2020). Such results show the importance of paying attention to these MNVs as highlighted by McLaren et al. (2016): “Current annotation tools, including the VEP, annotate each input variant independently, without considering the potential compound effects of combining alternate alleles across multiple variant loci.”

FIGURE 3

Example of an Erroneously Predicted Stop-Gained

As an example of erroneously predicted stop-gained, we present the case of the SLC27A4 gene, which is located on the reverse strand of chicken chromosome 17 (ENSGALG00000004965) (Figure 4A). In this gene, two SNPs rs316701182 and rs15031398, already reported in the Ensembl SNP database (Ensembl, 2018), were respectively, predicted as a stop-gained variant (TGA; stop-gained) and a synonymous variant (CGC; arginine) when compared to the reference haplotype (CGA; arginine) (Figure 4B). These SNPs were present in the FLLL population with frequencies > 20% and interestingly with contrasted frequencies between FL (fat line) and LL (lean line), two subpopulations divergently selected for adipose tissue weight (Leclercq et al., 1980). The rs15031398 SNP is absent in FL (Figure 4B); in the LL population in which we observed both SNPs (Figure 4C), we did not find any animal with the TGA (stop-gained) haplotype (composed of one variant only), with the rs316701182 T variant being always associated with the rs15031398 C variant within the “TGC” MNV. The absence of TGA (stop-gained) haplotype is consistent with several SLC27A4-knockout mouse studies which report prenatal lethality (Gimeno et al., 2003) or neonatal lethality (Herrmann et al., 2003; Moulson et al., 2003; Lin et al., 2010; Tao et al., 2012). The SLC27A4 gene codes fatty acid transport protein 4 (FATP4), which is particularly involved in the uptake of long-chain fatty acids (LCFAs); this gene is highly expressed in various chicken tissues as shown in Figure 4D with an expression > 10 TPM in the liver, ovary, optical system, skin, and intestine (ileum). Interestingly, FATP4 is thought to play a major role in dietary fatty acid uptake in intestinal epithelial cells (Hirsch et al., 1998) and in physiological uptake across cell membranes of LCFAs, which are key metabolites for energy generation and storage; it is viewed as a target to prevent or reverse obesity (Hirsch et al., 1998; Schaffer, 2002). FATP4 could be then related to the lean phenotype of the LL population for two reasons. First, the “TGC” (cysteine) MNV haplotype is reported as a severe change by the SIFT software package compared to the reference “CGA” (arginine) haplotype, suggesting a severe impact on the FATP4 protein function. Second, this “TGC” MNV haplotype is absent in FL birds, whereas it is frequent (42%) in 12 LL birds, with a higher frequency than expected (Figure 4C). We confirmed these results by extending this analysis to 58 birds (29 birds per line) using PCR amplification of the region of interest followed by Sanger sequencing. No rs15031398 was identified in the FL line. In the LL line, we observed 12 birds carrying the “TGC” MNV haplotype (three homozygous and nine heterozygous) and no bird with the TGA (stop-gained) haplotype. These results suggest a strong but not lethal impact of the MNV haplotype on the FATP4 protein function, which could then participate to the lean phenotype of the LL line. However, a genetic association study is needed to support a potential causal link between the FATP4 dysfunctional MNV and a low adiposity in the LL line compared to the FL line.

FIGURE 4

Jehl et al., 2020).

Conclusion

We have shown that MNVs represent an important class of genetic variations since they have a significant impact on polymorphism functional interpretation with roughly 40% of MNVs in our dataset inducing reannotation. These reannotations show a decreased impact severity of MNVs when compared to their constituent SNPs, at least for the stop-gained category. As previously demonstrated in human studies, our results in chicken demonstrate the value of haplotype-aware variant annotation and the interest to consider MNVs in coding region particularly when focusing on severe functional consequences such as stop-gained. We illustrated such a case with an erroneous stop-gained annotation found in the chicken SLC27A4 gene.

References

• 1

  BojanowskiM. (2020). mbojan/alluvial. R. Available online at: https://github.com/mbojan/alluvial. (accessed 8 Jan 2021).

  • Google Scholar
  • View reference in article

• 2

  ChengS.-J.ShiF.-Y.LiuH.DingY.JiangS.LiangN.et al (2017). Accurately annotate compound effects of genetic variants using a context-sensitive framework.Nucleic Acids Res.45:e82. 10.1093/nar/gkx041

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  CingolaniP.PlattsA.WangL. L.CoonM.NguyenT.WangL.et al (2012). A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3.Fly (Austin)680–92. 10.4161/fly.19695

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  DanecekP.McCarthyS. A. (2017). BCFtools/csq: haplotype-aware variant consequences.Bioinformatics332037–2039. 10.1093/bioinformatics/btx100

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  DePristoM. A.BanksE.PoplinR.GarimellaK. V.MaguireJ. R.HartlC.et al (2011). A framework for variation discovery and genotyping using next-generation DNA sequencing data.Nat. Genet.43491–498.

  • Google Scholar
  • View reference in article

• 6

  Ensembl. (2018). dbSNP - Gallus Gallus 5 - V94.ftp://ftp.ensembl.org/pub/release-94/variation/vcf/gallus_gallus(accessed December 11, 2020).^∗

  • Google Scholar
  • View reference in article

• 7

  GimenoR. E.HirschD. J.PunreddyS.SunY.OrtegonA. M.WuH.et al (2003). Targeted deletion of fatty acid transport protein-4 results in early embryonic lethality.J. Biol. Chem.27849512–49516. 10.1074/jbc.m309759200

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  HerrmannT.van der HoevenF.GröneH.-J.StewartA. F.LangbeinL.KaiserI.et al (2003). Mice with targeted disruption of the fatty acid transport protein 4 (Fatp 4, Slc27a4) gene show features of lethal restrictive dermopathy.J. Cell Biol.1611105–1115. 10.1083/jcb.200207080

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  HirschD.StahlA.LodishH. F. (1998). A family of fatty acid transporters conserved from mycobacterium to man.Proc Natl Acad Sci U S A.958625–8629. 10.1073/pnas.95.15.8625

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  JehlF.DegalezF.BernardM.LecerfF.LagoutteL.DésertC.et al (2021). RNA-Seq data for reliable SNP detection and genotype calling: Interest for coding variant characterization and cis-regulation analysis by allele-specific expression in livestock species. Front. Genet.12:655707. 10.3389/fgene.2021.655707

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  JehlF.MuretK.BernardM.BoutinM.LagoutteL.DésertC.et al (2020). An integrative atlas of chicken long non-coding genes and their annotations across 25 tissues.Sci. Rep.10:20457.

  • Google Scholar
  • View reference in article

• 12

  KhanW.Varma SaripellaG.LudwigT.CuppensT.ThibordF.GéninE.et al (2018). MACARON: a python framework to identify and re-annotate multi-base affected codons in whole genome/exome sequence data.Bioinformatics343396–3398. 10.1093/bioinformatics/bty382

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  LaiZ.MarkovetsA.AhdesmakiM.ChapmanB.HofmannO.McEwenR.et al (2016). VarDict: a novel and versatile variant caller for next-generation sequencing in cancer research.Nucleic Acids Res.44:e108. 10.1093/nar/gkw227

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  LeclercqB.BlumJ. C.BoyerJ. P. (1980). Selecting broilers for low or high abdominal fat: Initial observations.Br. Poult. Sci.21107–113. 10.1080/00071668008416644

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  LekM.KarczewskiK. J.MinikelE. V.SamochaK. E.BanksE.FennellT.et al (2016). Analysis of protein-coding genetic variation in 60,706 humans.Nature536285–291.

  • Google Scholar
  • View reference in article

• 16

  LinM.-H.ChangK.-W.LinS.-C.MinerJ. H. (2010). Epidermal hyperproliferation in mice lacking fatty acid transport protein 4 (FATP4) involves ectopic EGF receptor and STAT3 signaling.Dev. Biol.344707–719. 10.1016/j.ydbio.2010.05.503

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  McKennaA.HannaM.BanksE.SivachenkoA.CibulskisK.KernytskyA.et al (2020). The Genome Analysis Toolkit: A Mapreduce Framework For Analyzing Next-Generation Dna Sequencing Data. Available online at: https://genome.cshlp.org/content/20/9/1297.long. (accessed 18 Sep 2020).

  • Google Scholar
  • View reference in article

• 18

  McLarenW.GilL.HuntS. E.RiatH. S.RitchieG. R. S.ThormannA.et al (2016). The ensembl variant effect predictor.Genome Biol.17:122.

  • Google Scholar
  • View reference in article

• 19

  MoulsonC. L.MartinD. R.LugusJ. J.SchafferJ. E.LindA. C.MinerJ. H. (2003). Cloning of wrinkle-free, a previously uncharacterized mouse mutation, reveals crucial roles for fatty acid transport protein 4 in skin and hair development.Proc. Natl. Acad. Sci. U.S.A.1005274–5279. 10.1073/pnas.0431186100

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  NgP. C.HenikoffS. (2003). SIFT: predicting amino acid changes that affect protein function.Nucleic Acids Res.313812–3814. 10.1093/nar/gkg509

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  PagèsH.AboyounP.GentlemanR.DebRoyS. (2021). Biostrings: Efficient Manipulation Of Biological Strings. Bioconductor Version: Release (3.12).10.18129/B9.bioc.Biostrings

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  RobinsonJ. T.ThorvaldsdóttirH.WincklerW.GuttmanM.LanderE. S.GetzG.et al (2011). Integrative genomics viewer.Nat. Biotechnol.2924–26. 10.1038/nbt.1754

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  SchafferJ. E. (2002). Fatty acid transport: the roads taken.Am. J. Physiol. Endocrinol. Metab.282E239–E246.

  • Google Scholar
  • View reference in article

• 24

  SimN.-L.KumarP.HuJ.HenikoffS.SchneiderG.NgP. C. (2012). SIFT web server: predicting effects of amino acid substitutions on proteins.Nucleic Acids Res.40W452–W457.

  • Google Scholar
  • View reference in article

• 25

  SzklarczykD.GableA. L.LyonD.JungeA.WyderS.Huerta-CepasJ.et al (2019). STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets.Nucleic Acids Res.47D607–D613.

  • Google Scholar
  • View reference in article

• 26

  TaoJ.KosterM. I.HarrisonW.MoranJ. L.BeierD. R.RoopD. R.et al (2012). A spontaneous fatp4/scl27a4 splice site mutation in a new murine model for congenital ichthyosis.PLoS One.7:e50634. 10.1371/journal.pone.0050634

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  Van der AuweraG. A.CarneiroM. O.HartlC.PoplinR.Del AngelG.Levy-MoonshineA.et al (2013). From fastq data to high-confidence variant calls: the genome analysis toolkit best practices pipeline.Curr. Protoc. Bioinforma.4311.10.1–11.10.33.

  • Google Scholar

• 28

  WangK.LiM.HakonarsonH. (2010). ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data.Nucleic Acids Res.38:e164. 10.1093/nar/gkq603

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29

  WangQ.Pierce-HoffmanE.CummingsB. B.AlföldiJ.FrancioliL. C.GauthierL. D.et al (2020). Landscape of multi-nucleotide variants in 125,748 human exomes and 15,708 genomes.Nat. Commun.11:2539.

  • Google Scholar

• 30

  WeiL.LiuL. T.ConroyJ. R.HuQ.ConroyJ. M.MorrisonC. D.et al (2015). MAC: identifying and correcting annotation for multi-nucleotide variations.BMC Genomics.16:569. 10.1186/s12864-015-1779-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar