Published online 10 September 2019 Nucleic Acids Research, 2019, Vol. 47, No. 21 e136 doi: 10.1093/nar/gkz776 Assessment of kinship detection using RNA-seq data Natalia Blay1,2,3 , Eduard Casas1,2,4 , Iván Galván-Femenı́a1,5 , Jan Graffelman6,7 , Rafael de Cid1,5 and Tanya Vavouri 1,2,* 1 Program for Predictive and Personalized Medicine of Cancer, Germans Trias i Pujol Research Institute (PMPPC-IGTP), Badalona 08916, Spain, 2 Josep Carreras Leukaemia Research Institute (IJC), Campus ICO-Germans Trias i Pujol, Universitat Autònoma de Barcelona, Badalona 08916, Spain, 3 Masters Programme in Bioinformatics and Biostatistics, Universitat Oberta de Catalunya (UOC), Barcelona 08035, Spain, 4 Doctoral Programme in Biomedicine, Universitat de Barcelona, Barcelona 08007, Spain, 5 Genomes for Life - GCAT lab Group - Germans Trias i Pujol Research Institute, Can Ruti Campus, Ctra de Can Ruti, Camı́ de les Escoles s/n, Badalona, Barcelona 08916, Spain, 6 Department of Statistics and Operations Research Universitat Politècnica de Catalunya, Barcelona 08028, Spain and 7 Department of Biostatistics, University of Washington, Seattle, WA 98105-946, USA Received February 16, 2019; Revised August 23, 2019; Editorial Decision August 28, 2019; Accepted August 29, 2019 ABSTRACT (2), to understand the effects of disease-associated muta- tions (e.g. (3)), to study the mechanisms of epigenetic in- Analysis of RNA sequencing (RNA-seq) data from heritance of phenotypic traits (4), among others. In these related individuals is widely used in clinical and studies, multiple RNA samples are extracted from differ- molecular genetics studies. Prediction of kinship ent related or unrelated individuals and they are processed from RNA-seq data would be useful for confirming in parallel. Knowing whether samples are from related in- the expected relationships in family based studies dividuals and their exact relationship is important. For ex- and for highlighting samples from related individu- ample, in family studies, matching samples to individuals als in case-control or population based studies. Cur- in a pedigree allows the correct association of phenotypic rently, reconstruction of pedigrees is largely based traits to expression patterns. In case-control gene expres- on SNPs or microsatellites, obtained from genotyp- sion studies, the inclusion of samples from related individ- ing arrays, whole genome sequencing and whole ex- uals leads to a mis-specified covariance structure and there- fore an inflated type-1 error rate when testing for associa- ome sequencing. Potential problems with using RNA- tion. Knowing the true relatedness of the samples allows seq data for kinship detection are the low proportion researchers to remove unwanted related samples or use a of the genome that it covers, the highly skewed cov- method that models expression taking into account relat- erage of exons of different genes depending on ex- edness. In population studies, biased recruitment schemes pression level and allele-specific expression. In this can enrich datasets with cryptic relationships (5). This is a study we assess the use of RNA-seq data to detect common scenario in association analyses and it could also kinship between individuals, through pairwise iden- happen for RNA-seq collections from unrelated individu- tity by descent (IBD) estimates. First, we obtained als. Currently, researchers check the correctness of RNA- high quality SNPs after successive filters to minimize seq sample relatedness indirectly because it is so far un- the effects due to allelic imbalance as well as errors clear whether kinship can be detected directly from RNA- in sequencing, mapping and genotyping. Then, we seq data. During sample and data processing, sample mislabeling used these SNPs to calculate pairwise IBD estimates. can occur due to human error. This is detrimental for down- By analysing both real and simulated RNA-seq data stream analysis, especially for family studies, and has been we show that it is possible to identify up to second estimated to affect at least 4% of published samples (6,7). degree relationships using RNA-seq data of even low Given the importance and pervasiveness of this problem, to moderate sequencing depth. there are numerous programs and methodologies that deal with sample mislabeling in sequencing data by compar- ing paired samples from the same individual (8–11), com- INTRODUCTION paring the annotated sex with expression of sex-specific RNA sequencing is used in genetics for a variety of pur- genes (12) or using a heuristic data perturbation strategy poses such as to test the heritability of gene expression traits (7). From genotyping, whole genome sequencing (WGS) (1), to search for mutations causing mendelian disorders and whole exome sequencing (WES) datasets, there are nu- * To whom correspondence should be addressed. Tel: +34 935543078; Fax: +34 934651472; Email: tvavouri@carrerasresearch.org  C The Author(s) 2019. Published by Oxford University Press on behalf of Nucleic Acids Research. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. e136 Nucleic Acids Research, 2019, Vol. 47, No. 21 PAGE 2 OF 9 merous methods that can predict kinship from single nu- two unrelated individuals (27) as well as targeted transcript cleotide polymorphisms (SNPs) and confirm that samples sequencing (Ion AmpliSeq, Life Technologies) from 7 un- correspond to the labelled individuals in a pedigree or that related individuals and a pair of siblings (28). Second, we they are genetically unrelated. However, none of these meth- retrieved genetic variation data from six pairs of first and ods sufficiently address labeling mix-ups or unreported re- second degree relatives from the CDX population of the latedness of samples when only RNA-seq data is available, 1000 Genomes Project (29). From this data we then sim- for each individual. Furthermore, for researchers interested ulated RNA-seq reads (the method of simulation of RNA- in analyzing only expression data, current methods require seq reads is described below). Third, we retrieved genetic them to retrieve and process in parallel at least one addi- variation data from unrelated individuals from the 1000 tional data file (WGS, WES or genotyping array data) per Genomes Project (29) from which we simulated families (the individual, so as to confirm that RNA-seq and DNA se- method of simulation of genotypes of family members is de- quencing data contain the same nucleotide variants and scribed below) and RNA-seq reads. then use the DNA sequencing data to extract genotypes from which to confirm kinship. Even when available, the re- Simulation of genotypes of family members trieval and analysis of WGS or WES just to check for errors in RNA-seq sample relatedness is a cumbersome and often To simulate the genotypes of family members we used SNP very time consuming task. data from the 1000 Genomes Project (29). We selected unre- There are several methods to predict or confirm familial lated individuals as founders and used their haplotype data relationships and build heritability estimates from genetic to simulate the rest of the family. We simulated data ac- data. Although these methods were initially based on mi- cording to different pedigrees in order to assess the effect crosatellites (13) and amplified fragment length polymor- of pedigree complexity. Family simulations were carried out phism (AFLPs) (14), nowadays they usually rely on SNPs. following Mendelian laws and taking into account linkage Among the different methodologies for kinship detec- disequilibrium (assuming 1 cM per Mb) with a custom R tion, we can find methods that exclude impossible parent- script. offspring relationships (15–17), methods that calculate kin- ship coefficients and/or identical-by-descent (IBD) proba- Simulation of RNA-seq reads bilities (18,19) and likelihood methods (20,21). RNA-seq To simulate gene expression data from real or simulated data can be used to identify SNPs (22). Yet, the use of RNA- individuals, we generated paired-end RNA-seq reads with seq data for kinship detection and pedigree reconstruction flux-simulator v1.2.1 (30). We used two genome fasta files has not been properly assessed so far. Some of the concerns per individual (one per haplotype, including the SNPs to with using this type of data for kinship detection are the the reference genome with GATK v3.8 FastaAlternateRef- limited number of SNPs adequately covered by reads and erenceMaker (31)) and the expression profile of B lympho- allelic imbalance (23,24) that may hinder genotype calling cytes (custom .pro file with RPKMx150 obtained from the by masking one of the two alleles at many variant positions founders of CEPH/UTAH family 1463 (26) to get a se- along the genome. quencing depth of 40M reads per individual). We ran the Here, we report the results of our assessment of the utility simulation with library preparation and sequencing steps (- of human RNA-seq data for kinship detection and pedigree ls options) for the first haplotype. For the second haplotype reconstruction. Our input data is a set of raw RNA-seq data we used the library (.lib file) of the first haplotype and ran files and a set of known common variant positions. First, we only the sequencing step (-s option) to make sure the same filter both the RNA-seq data and the known variants to ob- genes are expressed in both haplotypes. tain reliable predicted genotypes at common variant posi- tions for each RNA sample. We then use the predicted geno- Data filtering, read mapping and kinship detection types to detect and represent familial relationships by esti- mating pairwise probabilities of identity-by-descent (IBD) We aligned empirical and simulated RNA-seq reads to the (18,25). By analysing human empirical as well as simulated reference genome hg19 using HISAT2 v2.1.0 (32). For map- RNA-seq data we show that the estimated IBD probabil- ping, we used the HISAT2 index for the human reference ities allow kinship detection and pedigree reconstruction, genome plus transcripts retrieved from ftp://ftp.ccb.jhu.edu/ detecting up to second degree relationships even with low pub/infphilo/hisat2/data/grch37 tran.tar.gz (accessed: Nov sequencing depth. We propose kinship detection directly 8, 2018). We then used SAMtools to remove duplicates from RNA-seq data as a simple and direct method to de- (command markdup). We retrieved 14.8 million (M) com- tect sample relatedness. mon SNPs from the UCSC Genome Browser (dbSNP Build ID 150) (33). We used BEDtools intersect v2.26.0 (34) to MATERIALS AND METHODS identify and remove 38,572 SNPs in 83 imprinted genes (http://www.geneimprint.org/ accessed: 8 March 2018) and Retrieval of empirical data 8.2 M SNPs in repeats (RepeatMasker annotation down- We used three types of empirical data (Figure 1, Supplemen- loaded from the UCSC Genome Browser), ending up with tary Table S1). First, we used previously published RNA- 6.2 M SNPs. Genotypes were obtained with SAMtools seq data from four different studies that include transcrip- mpileup v1.7 (options -A -q 4 -t AD,DP) (35) and BCFtools tome sequencing from a 17-member, three-generation fam- call v1.4 (options -m - -O b -f GQ) (36), using uniquely ily (26), transcriptome sequencing from a trio (3), transcrip- mapping reads only. We used VCFtools v0.1.14 (36) to se- tome sequencing from a pair of first degree relatives and lect only those SNPs with a depth ≥ 10 (option –minDP PAGE 3 OF 9 Nucleic Acids Research, 2019, Vol. 47, No. 21 e136 Figure 1. Overview of data workflow for kinship detection and pedigree reconstruction using RNA-seq data. GQ: genotype quality, DP:depth, IBD: identity by descent, MAF: minor allele frequency. 10) and a genotype quality ≥ 20 (option –minGQ 20). To A B assess the accuracy of genotype prediction from RNA-seq data, we retrieved high confidence variant calls for two indi- viduals of CEPH/UTAH family 1463 (the two offspring of the four founders in the pedigree shown in Figure 2A) (37). The effect of the different SNP filtering steps on the num- ber of variants considered for IBD estimation, the accu- racy of genotype prediction (for individuals NA12877 and C NA12878 of CEPH/UTAH family 1463), kinship detection and pedigree reconstruction is shown in Supplementary Ta- ble S2. To determine pairwise IBD estimates we used PLINK v1.9 (38), considering only autosomal SNPs with a minor allele frequency ≥0.3 (–maf 0.3 option) to obtain an opti- mal separation between groups. Use of a higher minor al- lele frequency threshold led to high deviation of Z0, Z1 and Z2 from the theoretical values (Supplementary Figure S1). Representation of pairwise IBD estimates was done in R Figure 2. Structure of datasets used for the assessment of kinship detection using RNA-seq data. (A) Structure of the extended CEPH/UTAH family v3.5.0 (39) with the method described by Galván-Femenı́a 1463 with empirical RNA-seq data (26). Additional empirical RNA-seq et al. (40). Using the IBD estimates obtained in the previous datasets are included in Supplementary Figure S2. (B) List of real related step, we used PRIMUS v1.9.0 (41) to predict pairwise rela- pairs of individuals with simulated RNA-seq data. (C) Simulated fami- tionships and to reconstruct the whole pedigree, consider- lies with simulated RNA-seq data (pedigree types 1–4). Real individuals with simulated RNA-seq data are highlighted in grey. Simulated individu- ing as unrelated those pairs of individuals with a coefficient als with simulated RNA-seq data are highlighted in black. of relatedness lower than 0.2 (option -t 0.2). Sex data was also provided. The sex of the individuals from the real data was inferred from counting reads mapping on the Y chro- mosome and compared against the reported sex. e136 Nucleic Acids Research, 2019, Vol. 47, No. 21 PAGE 4 OF 9 RESULTS that we have analysed. These high confidence variant calls agreed with those predicted from RNA-seq data using our Analysis workflow for kinship detection using RNA-seq from pipeline at 99.3% and 98.8% of the common variant posi- related individuals tions that passed all filters. To assess how feasible it is to detect kinship using RNA-seq We used ternary diagrams to visualize the estimated data we used three types of data (Figures 1 and 2). First, we probabilities of sharing 0, 1 and 2 IBD alleles (referred here used six previously published sets of RNA-seq data from as Z0, Z1 and Z2 respectively) meaning that none, only one related and unrelated individuals from four different stud- or both alleles of the pair are inherited from the same recent ies that include samples from an extended three-generation common ancestor. In ternary diagrams, the coordinates of human family (26) (Figure 2A). Second, we used available each data point are the three IBD estimates for each pair SNP data from six pairs of first and second degree relatives of individuals in the dataset. We observed that the differ- from the CDX population of the 1000 Genomes Project ent relationships were grouped around their theoretical val- (29) from which we simulated RNA-seq data (Figure 2B). ues with no overlap between them. As expected, parent- Third, we used available SNP data from unrelated individ- offspring pairs are clustered at the Z1 corner (meaning that uals from the 1000 Genomes Project (29) from which we they share 1 IBD allele at all SNPs), unrelated individuals simulated the genotypes of family members and the corre- at the Z0 corner (they do not share any IBD alleles), second sponding RNA-seq samples (Figure 2C). We simulated four degree relationships on the Z0-Z1 axis (they share 1 IBD al- types of pedigrees to assess the effect of pedigree complexity. lele at half of their SNPs) and full-siblings in the middle of After data retrieval or simulation, we used the same analysis the graph (they can share from 0 to 2 IBD alleles per SNP) workflow which in brief consisted of mapping and filtering (Figure 3A, Supplementary Table S3). RNA-seq reads, calling and filtering genotypes of a set of The actual pedigree was also correctly reconstructed with filtered SNPs, calculating IBD estimates and finally visual- PRIMUS (41) from the estimated IBD values. In addition izing the IBD estimates in ternary diagrams or providing to the top scoring pedigree, which in this case corresponds the IBD estimates to the software PRIMUS for prediction to the correct one, PRIMUS also reports the predicted pair- of pairwise relationships and pedigree reconstruction. wise relationships. In this case it did not correctly identify some of the pairwise relationships; the mis-identified pairs were 2% of full sibling pairs (one pair was classified as sec- Kinship detection using empirical human RNA-seq data ond degree relatives), 32% of second degree pairs and 9% of We mapped and filtered the raw RNA-seq data from B- unrelated pairs of individuals. lymphocytes from members of the CEPH/UTAH family In addition to the data from the 17-member family, we 1463 (26) and in parallel filtered the initial set of known assessed kinship detection for a family trio (mother, father, common human variants as described in the methods. daughter) using transcriptome sequencing data from car- Reads were mapped to the reference human genome. This diomyocytes differentiated from induced pluripotent stem step excluded on average 12% of the initial reads (range cells (3), for a parent-offspring pair and unrelated indi- 6–24%, depending on the individual) that could not be viduals using transcriptome sequencing data from whole mapped. We then removed PCR duplicates as they can po- blood and peripheral blood mononuclear cells (27) and for tentially amplify sequencing errors. This step led to the re- a pair of siblings and seven unrelated individuals using tar- moval of an average of 21% of the reads (range 17–30%). geted RNA sequencing from naive CD4+ T-cells and LPS- We also removed multi-mapping reads to avoid mapping stimulated monocytes (28) (Supplementary Figure S2A–E). errors, excluding on average 4.6% of the reads (range 4.1– For whole transcriptome sequencing, the number of SNPs 5.3%). We filtered the known common human variants; we used for pairwise comparisons and IBD estimation ranged removed SNPs in repeats and in imprinted genes (55% and from 1,311 to 5,876 SNPs, depending on the pair of in- 0.26% of SNPs excluded, respectively) that could bias geno- dividuals analysed. For targeted RNA-seq samples, gener- typing. Having obtained a filtered set of RNA-seq reads (a ated using the Ion AmpliSeq Transcriptome Human Gene mean of 34M reads per individual) and SNPs (6.2 M SNPs), Expression Kit (Life Technologies) that targets on aver- we proceeded to call and filter genotypes, obtaining 633,636 age 150 bp per human gene, the number of SNPs cov- SNPs that were covered by reads in at least one individual. ered and passing all filters were significantly fewer, as ex- To exclude those SNPs with low coverage we selected a min- pected. Specifically, from Ion AmpliSeq data, the number imum depth of 10 reads per SNP, and to avoid genotyping of SNPs used for pairwise comparisons and IBD estima- errors we filtered out those SNPs with a genotype quality tion ranged between 95 and 193. Samples from related indi- lower than 20, which means an error rate of 0.01 (in this step viduals (parent-offspring pairs and siblings) are clearly sep- we excluded 76% of the SNPs). Finally, we considered those arated from replicates (tightly clustering at the Z2 corner) SNPs with a minor allele frequency of at least 0.3 (exclud- and from unrelated individuals in all datasets, including in ing 89.5% of the SNPs). We ended the filtering process with the targeted RNA-seq datasets (Supplementary Figure S2). 16,004 SNPs present in at least one individual. The num- In these datasets, although samples from unrelated individ- ber of SNPs shared between pairs of individuals that were uals are clearly separated from samples from first-degree used for downstream analyses was between 5,657 and 8,600, relatives, their IBD estimates are further from their theo- depending on the number of missing genotypes observed retical values, something that was improved by removing for each pair of individuals. High confidence variant calls replicates (data not shown). Relatedness was correctly in- are available for two individuals (individuals with identifiers ferred by PRIMUS, however for the transcriptome sequenc- NA12877 and NA12878) of the CEPH/UTAH family 1463 ing data from whole blood and the targeted RNA-seq data PAGE 5 OF 9 Nucleic Acids Research, 2019, Vol. 47, No. 21 e136 A B C Figure 3. Ternary diagrams of IBD estimates for (A) CEPH/UTAH family 1463, (B) simulated RNA-seq data from real related pairs of individuals and (C) simulated RNA-seq data from simulated related individuals from different pedigree types (type 1–type 4). Note that in ternary diagrams data points that have very similar Z0, Z1 and Z2 values overlap and may not be visible. Ternary diagrams for five additional empirical RNA-seq datasets are shown in Supplementary Figure S2. from peripheral blood mononuclear cells other high scoring subsampled 37M reads (the number of reads for the individ- pedigrees were also suggested. ual with the lowest coverage) to homogenize the sequencing According to these results, we can conclude that kinship depth for all individuals and then also 30M, 20M, 10M, 8M, detection is possible with RNA-seq data, allowing us to de- 6M, 4M and 2M reads from each sample to determine the tect and represent first and second degree relationships. Al- minimum sequencing depth for kinship detection and pedi- though the sensitivity decreases as the degree of relatedness gree reconstruction. increases, pedigree reconstruction is robust. As an indica- Visual inspection of estimated probabilities using ternary tion of the time required to perform kinship detection, the diagrams (Figure 4) revealed that there was no overlap full RNA-seq analysis of nine samples from the trio (fa- between second degree relatives and unrelated individuals ther, mother and offspring with three replicates per individ- when using RNA-seq data with 52M, 37M, 30M, 20M, ual from (3)) took 110 min using an Intel(R) Xeon(R) W- 10M and 6M reads. In the simulation of 8M reads we ob- 2145 CPU 3.70GHz workstation with eight cores and 64GB served a slight overlap between second degree relatives and of memory (parallel processing was used only for mapping unrelated individuals. In the simulation of 4M reads, there reads). This includes the time needed to map the raw reads was overlap between second degree and unrelated pairs of (∼45 min). individuals and also between second degree and full sib- lings. With 2M reads all groups overlapped except parent- Impact of sequencing depth on kinship detection offspring and unrelated pairs of individuals. Pedigree reconstruction was only possible for RNA-seq The sequencing depth of the real RNA-seq data from the ex- data with 52M, 37M, 30M and 20M reads because, when tended human family ranged from 37M to 67M reads, with decreasing the sequencing depth, the probability to share a mean depth of 52M. To understand the effect of sequenc- one IBD allele (Z1) of some full siblings decreased way be- ing depth on kinship detection, we randomly subsampled low the expected values (e.g. Z0 = 0.6741, Z1 = 0.0712, Z2 different numbers of reads from the FASTQ input files. We = 0.2547). e136 Nucleic Acids Research, 2019, Vol. 47, No. 21 PAGE 6 OF 9 of first and second degree relatives from the CDX popu- lation of the 1000 Genomes Project (Figure 2B). We simu- lated paired-end RNA-seq reads with the expression profile of B-lymphocytes. We obtained this gene expression pro- file from the mean of the number of reads mapped to a gene per gene length in kilobases per million mapped reads (RPKM) of the founders of the CEPH/UTAH family 1463. Then, we multiplied those values by a factor of 150 to sim- ulate the desired expression profile with a sequencing depth of 40M reads. As expected, the simulated RNA-seq data is highly correlated with the real data (Supplementary Figure S3)––Pearson’s correlation estimates are on average 0.918 (range 0.913–0.925). We used the same workflow described for the empirical RNA-seq data. The resulting number of SNPs used for the pairwise comparisons was between 4,501 and 4,632. In agreement with the results obtained from the real RNA-seq data, visual inspection of estimated probabili- ties in ternary diagrams revealed that all relationships are clearly separated (Figure 3B). Similarly, the pedigree recon- struction program was able to correctly classify and repre- sent all relationships in a pedigree. We conclude that kin- ship detection using RNA-seq data is possible using differ- ent datasets; for entire families as well as for pairs of related individuals. Kinship detection using simulated genotypes of family mem- bers with simulated RNA-seq data We then asked whether the complexity of the pedigree could affect kinship detection and pedigree reconstruction. To ad- dress this, we generated additional RNA-seq data from dif- ferent types of pedigrees, this time simulating both gene ex- pression and also the genotypes of individuals in the pedi- Figure 4. Ternary diagrams of IBD estimates for each sequencing depth grees. We used genotypes of real unrelated individuals from (top) and number of SNPs used for pairwise comparison per sequencing the 1000 Genomes Project as founders (29) (Figure 2C, depth (bottom). Green color indicates that the pedigree was correctly re- constructed, yellow colour indicates that there was no overlapping between genotypes of real individuals shown in gray, simulated geno- groups but the pedigree reconstruction was not possible, and red color in- types of offspring shown in black). We simulated families dicates that there was an overlapping between groups. The box with a se- using data from eight different populations from the 1000 quencing depth of 52 is the one for the original data, being the mean for Genomes Project: GBR, KHV, IBS, LWK, CLM, CDX, the different sequencing depths (from 37M to 67M reads). PEL and ACB; and four pedigree structures including dif- ferent degrees of relationship (Figure 2C): Taken together, these results demonstrate that sequenc- Type 1: A 16-member pedigree with the same relationships ing depth affects the number of SNPs used for pairwise as CEPH/UTAH family 1463. comparison, but overall its effect on kinship detection is Type 2: A 16-member pedigree with all first and second de- small, providing acceptable IBD estimates (no overlapping gree relationships from Supplementary Table S3. between relationships) even with 6M reads, where the num- Type 3: A 16-member pedigree with first, second and third ber of SNPs is as low as 811–1,282 (Figure 4). As IBD es- degree relationships. timates with 8M reads were not as good, we conclude that Type 4: A 12-member pedigree with all first, second and the minimum sequencing depth to be able to detect kinship third degree relationships from Supplementary Table S3. in RNA-seq data through IBD probabilities is 10M reads. In contrast, the minimum sequencing depth to be able to We made 32 simulations, one for each pedigree type and reconstruct the correct pedigree is 20M reads. population. Sex was not simulated as sex chromosomes are not used for IBD estimation. After data filtering, we ob- served that the number of SNPs used for pairwise compar- Kinship detection using empirical genotypes of related indi- ison was higher in African populations (ACB and LWK) viduals with simulated RNA-seq data ranging from 5,180 to 5,974 and lower in Asian popula- To further assess the feasibility of estimating relationships tions (KHV and CDX) ranging from 4,289 to 5,084 (Figure from RNA-seq samples from related individuals, we sim- 5B). To further assess the accuracy of our genotype predic- ulated RNA-seq data from real related individuals with tion pipeline, we compared the genotypes of individuals in a known genotypes. We used variation data from six pairs simulated pedigree (type 4 pedigree with founders from IBS PAGE 7 OF 9 Nucleic Acids Research, 2019, Vol. 47, No. 21 e136 A B Figure 5. Pedigree reconstruction using RNA-seq data. (A) The matrices show the correctly reconstructed simulated pedigrees with our default workflow (left panel), with five additional unrelated individuals (middle panel) and with a higher relatedness cutoff (0.375 instead of 0.2) (right panel). (B) The boxplots summarize the distributions of the number of SNPs used for pairwise comparisons and IBD estimation for each pair of individuals in each simulated pedigree. population) to the predicted genotypes and found them to four families, reconstructing correctly all type 1, seven type agree at 98.3–99.2% of the common variant sites that passed 2, all type 4 pedigrees but none of type 3 (Figure 5A, Mid- all filters. dle Panel). We conclude that adding data from unrelated Visual inspection of IBD probabilities (Figure 3C) re- individuals would be useful in those cases where the recon- vealed that first degree relationships were clearly separated struction is not possible. from other relationships, second degree relationships were We asked whether it was possible to correctly predict the separated from unrelated individuals in all simple pedigrees pedigree structure using only the most related individuals (type 1 simulated pedigrees) and some type 2, 3 and 4 simu- in complex pedigrees by using a stricter relatedness cutoff. lated pedigrees. In contrast, third degree relationships over- We reran PRIMUS with a higher relatedness cutoff (option lapped unrelated pairs of individuals in all pedigrees where -t 0.375, instead of the previously used -t 0.2) and found they were present (Figure 3C, pedigree type 3–type 4). that type 1 and type 2 pedigrees were reconstructed equally Reconstruction of the actual pedigree was possible for well but also that seven out of eight type 3 pedigrees and all all simulated type 1 pedigrees, seven type 2 pedigrees and type 4 pedigrees could also be correctly reconstructed (Fig- four type 4 pedigrees (one of them was not the top scoring ure 5A, right panel). We conclude that a more restrictive one, but the fourth in the PRIMUS output), but for none relatedness cutoff works better than adding unrelated indi- of the type 3 pedigrees (Figure 5A, left panel). Although viduals in those cases where most or all of the individuals of PRIMUS was not able to reconstruct all the simulated pedi- the pedigree are connected through a first degree relative. grees, it was able to identify all first degree relatives, 35% of Last, we tested whether IBD estimation and kinship de- second degree relatives (in type 1 pedigrees this percentage tection could be used to identify RNA-seq samples from re- was 68%), 12% of third degree relatives, and 93% of unre- lated individuals in a dataset consisting mostly of samples lated individuals (Supplementary Figure S4). from unrelated individuals. To test this we retrieved geno- We conclude that family structure affects kinship detec- types for twenty unrelated individuals of the IBS popula- tion, providing good IBD estimates for simple families with tion from the 1000 Genomes Project (29) (Supplementary up to second degree relatives, but in pedigrees with third Table S1) and we used them to simulate RNA-seq data as degree relatives there is overlap between different relation- done previously. We then generated five different RNA-seq ships. Accordingly, entire pedigree reconstruction works datasets, each consisting of the twenty simulated RNA-seq better for simple families, getting worse when the number samples from unrelated individuals from the 1000 Genomes of third degree related pairs increases (type 3 pedigree has Project (29) and a pair of the simulated RNA-seq samples 25 pairs of third degree related individuals while type 4 pedi- from the previously generated type 4 pedigree (with IBS gree has six). Population specific differences in genetic vari- genotypes as founders) (Supplementary Figure S5). The es- ability lead to different number of SNPs that pass all filters, timated IBD probabilities clearly separated a pair of RNA- however these differences do not affect the results. seq samples from second degree relatives (grandmother- We wondered whether the addition of unrelated individ- grandson) from samples from twenty unrelated individ- uals would improve kinship detection. Presumably, they in- uals from the same population. Samples from third de- clude more variability and make related individuals look gree relatives (great-grandmother and great-grandson and more similar to each other. To test this, we added five un- first cousins) had estimated IBD probabilities that were related individuals from the same population to each simu- marginally separated from those of unrelated individuals. lated family. The result was just as good for all families that Samples from fourth degree relatives were indistinguishable were previously correctly reconstructed but it improved for from those from unrelated individuals. We tested whether e136 Nucleic Acids Research, 2019, Vol. 47, No. 21 PAGE 8 OF 9 more stringent filtering of variants by removing variants ACKNOWLEDGEMENTS falling in segmental duplications would improve the sepa- We thank the IGTP High Performance Computing Plat- ration of third and fourth degree relatives from unrelated form (IGTP HPC) for computing resources and Iñaki M individuals but this was not the case (Supplementary Fig- de Ilarduya and Lloyd Goodman for technical support. ure S6). We conclude that IBD probabilities estimated from RNA-seq data can be used to confidently identify samples from up to second degree relatives and possibly flag samples FUNDING from third degree relatives within larger RNA-seq datasets Spanish Ministry of Economy and Competitive- from unrelated individuals. ness [BFU2015-70581 and ADE 10/00026]; Catalan Agency for Management of University and Research Grants [2017 SGR 1262 and 2017 SGR 529]; CERCA DISCUSSION Programme/Generalitat de Catalunya; Research at the IJC We have shown here that kinship detection based on es- was supported by the ‘La Caixa’ Foundation, the Fundació timates of identity by descent probabilities using RNA- Internacional Josep Carreras and Celgene Spain. Funding seq data is possible allowing the detection of up to second for open access charge: Catalan Agency for Management degree relationships. In addition, we have shown that the of University and Research Grants. actual pedigrees can be fully reconstructed and, although Conflict of interest statement. None declared. pedigree reconstruction works better for simple pedigrees, some pedigrees with third degree relatives can also be recon- REFERENCES structed. Furthermore, using simulations, we have shown 1. Price,A.L., Helgason,A., Thorleifsson,G., McCarroll,S.A., Kong,A. that the ability to detect kinship between individuals is not and Stefansson,K. (2011) Single-tissue and cross-tissue heritability of limited by RNA-seq sequencing depth since it is still pos- gene expression via identity-by-descent in related or unrelated sible at a depth significantly below what is considered ac- individuals. PLoS Genet., 7, e1001317. 2. Gonorazky,H.D., Naumenko,S., Ramani,A.K., Nelakuditi,V., ceptable for gene expression analyses (usually 50M reads). Mashouri,P., Wang,P., Kao,D., Ohri,K., Viththiyapaskaran,S., Furthermore, we have shown that samples from first de- Tarnopolsky,M.A. et al. (2019) Expanding the boundaries of RNA gree relatives can be distinguished from replicates and sam- sequencing as a diagnostic tool for rare mendelian disease. Am. J. ples from unrelated individuals even with targeted RNA-seq Hum. Genet., 104, 466–483. 3. Gifford,C.A., Ranade,S.S., Samarakoon,R., Salunga,H.T., de data. Reconstruction of full pedigrees requires higher se- Soysa,T.Y., Huang,Y., Zhou,P., Elfenbein,A., Wyman,S.K., Bui,Y.K. quencing depth but is still possible with at least 20M reads. et al. (2019) Oligogenic inheritance of a human heart disease An issue that remains to be investigated in future studies involving a genetic modifier. Science, 364, 865–870. is whether gene expression data from other species would 4. Heard,E. and Martienssen,R.A. (2014) Transgenerational epigenetic provide similar results to the ones from humans. Genetic inheritance: myths and mechanisms. Cell, 157, 95–109. 5. Voight,B.F. and Pritchard,J.K. (2005) Confounding from cryptic variability among different individuals from a particular relatedness in case-control association studies. PLoS Genet., 1, e32. species or a particular population could play an important 6. Broman,K.W., Keller,M.P., Broman,A.T., Kendziorski,C., role in kinship detection with higher genetic diversity be- Yandell,B.S., Sen,Ś. and Attie,A.D. (2015) Identification and tween individuals leading to better results. correction of sample Mix-Ups in expression genetic data: A case study. G3: Genes|Genomes|Genetics, 5, 2177–2186. In conclusion, we have shown here that RNA-seq data 7. Zych,K., Snoek,B.L., Elvin,M., Rodriguez,M., Van der Velde,K.J., can be used to identify samples from closely related indi- Arends,D., Westra,H.-J., Swertz,M.A., Poulin,G., Kammenga,J.E. viduals using estimated identity by descent probabilities cal- et al. (2017) reGenotyper: Detecting mislabeled samples in genetic culated from predicted genotypes at common variant posi- data. PLoS One, 12, e0171324. tions. It turns out that the effect of calling genotypes cor- 8. Lee,S., Lee,S., Ouellette,S., Park,W.-Y., Lee,E.A. and Park,P.J. (2017) NGSCheckMate: software for validating sample identity in rectly at the single nucleotide level due to allele-specific ex- next-generation sequencing studies within and across data types. pression is not large enough to impede kinship detection Nucleic Acids Res., 45, e103. from IBD estimates. We therefore recommend the estima- 9. Fort,A., Panousis,N.I., Garieri,M., Antonarakis,S.E., tion of IBD probabilities and visualization of the cluster- Lappalainen,T., Dermitzakis,E.T. and Delaneau,O. (2017) MBV: a method to solve sample mislabeling and detect technical bias in large ing of samples in ternary graphs or the use of pedigree re- combined genotype and sequencing assay datasets. Bioinformatics, construction programs such as PRIMUS as a quality con- 33, 1895–1897. trol step in studies that generate multiple RNA-seq samples 10. Huang,J., Chen,J., Lathrop,M. and Liang,L. (2013) A tool for RNA from related individuals and population based studies. sequencing sample identity check. Bioinformatics, 29, 1463–1464. 11. Yoo,S., Huang,T., Campbell,J.D., Lee,E., Tu,Z., Geraci,M.W., Powell,C.A., Schadt,E.E., Spira,A. and Zhu,J. (2014) MODMatcher: multi-omics data matcher for integrative genomic analysis. PLoS DATA AVAILABILITY Comput. Biol., 10, e1003790. 12. Toker,L., Feng,M. and Pavlidis,P. (2016) Whose sample is it anyway? The scripts used for kinship detection and pedigree recon- Widespread misannotation of samples in transcriptomics studies. struction are available through the GitHub repository http: F1000Res., 5, 2103. //github.com/vavouri-lab/RNA-kinship 13. Jones,A.G., Small,C.M., Paczolt,K.A. and Ratterman,N.L. (2010) A practical guide to methods of parentage analysis. Mol. Ecol. Resour., 10, 6–30. 14. Savelkoul,P.H., Aarts,H.J., de Haas,J., Dijkshoorn,L., Duim,B., SUPPLEMENTARY DATA Otsen,M., Rademaker,J.L., Schouls,L. and Lenstra,J.A. (1999) Amplified-fragment length polymorphism analysis: the state of an Supplementary Data are available at NAR Online. art. J. Clin. Microbiol., 37, 3083–3091. PAGE 9 OF 9 Nucleic Acids Research, 2019, Vol. 47, No. 21 e136 15. Thompson,E.A. and Meagher,T.R. (1987) Parental and sib 29. The 1000 Genomes Project Consortium (2015) A global reference for likelihoods in genealogy reconstruction. Biometrics, 43, 585–600. human genetic variation. Nature, 526, 68–74. 16. Calus,M.P.L., Mulder,H.A. and Bastiaansen,J.W.M. (2011) 30. Griebel,T., Zacher,B., Ribeca,P., Raineri,E., Lacroix,V., Guigó,R. Identification of Mendelian inconsistencies between SNP and and Sammeth,M. (2012) Modelling and simulating generic RNA-Seq pedigree information of sibs. Genet. Sel. Evol., 43, 34. experiments with the flux simulator. Nucleic Acids Res., 40, 17. Hayes,B.J. (2011) Efficient parentage assignment and pedigree 10073–10083. reconstruction with dense single nucleotide polymorphism data. J. 31. McKenna,A., Hanna,M., Banks,E., Sivachenko,A., Cibulskis,K., Dairy Sci., 94, 2114–2117. Kernytsky,A., Garimella,K., Altshuler,D., Gabriel,S., Daly,M. et al. 18. Thompson,E.A. (1975) The estimation of pairwise relationships. Ann. (2010) The Genome Analysis Toolkit: a MapReduce framework for Hum. Genet., 39, 173–188. analyzing next-generation DNA sequencing data. Genome Res., 20, 19. Weir,B.S., Anderson,A.D. and Hepler,A.B. (2006) Genetic 1297–1303. relatedness analysis: modern data and new challenges. Nat. Rev. 32. Kim,D., Langmead,B. and Salzberg,S.L. (2015) HISAT: a fast spliced Genet., 7, 771–780. aligner with low memory requirements. Nat. Methods, 12, 357–360. 20. Meagher,T.R. and Thompson,E. (1986) The relationship between 33. Sherry,S.T. (2001) dbSNP: the NCBI database of genetic variation. single parent and parent pair genetic likelihoods in genealogy Nucleic Acids Res., 29, 308–311. reconstruction. Theor. Popul. Biol., 29, 87–106. 34. Quinlan,A.R. and Hall,I.M. (2010) BEDTools: a flexible suite of 21. Hill,W.G., Salisbury,B.A. and Webb,A.J. (2008) Parentage utilities for comparing genomic features. Bioinformatics, 26, 841–842. identification using single nucleotide polymorphism genotypes: 35. Li,H., Handsaker,B., Wysoker,A., Fennell,T., Ruan,J., Homer,N., application to product tracing. J. Anim. Sci., 86, 2508–2517. Marth,G., Abecasis,G., Durbin,R. and 1000 Genome Project Data 22. Piskol,R., Ramaswami,G. and Li,J.B. (2013) Reliable identification of Processing Subgroup (2009) The Sequence Alignment/Map format genomic variants from RNA-seq data. Am. J. Hum. Genet., 93, and SAMtools. Bioinformatics, 25, 2078–2079. 641–651. 36. Danecek,P., Auton,A., Abecasis,G., Albers,C.A., Banks,E., 23. Deng,Q., Ramsköld,D., Reinius,B. and Sandberg,R. (2014) DePristo,M.A., Handsaker,R.E., Lunter,G., Marth,G.T., Sherry,S.T. Single-cell RNA-seq reveals dynamic, random monoallelic gene et al. (2011) The variant call format and VCFtools. Bioinformatics, expression in mammalian cells. Science, 343, 193–196. 27, 2156–2158. 24. Vigneau,S., Vinogradova,S., Savova,V. and Gimelbrant,A. (2018) 37. Eberle,M.A., Fritzilas,E., Krusche,P., Källberg,M., Moore,B.L., High prevalence of clonal monoallelic expression. Nat. Genet., 50, Bekritsky,M.A., Iqbal,Z., Chuang,H.-Y., Humphray,S.J., 1198–1199. Halpern,A.L. et al. (2017) A reference data set of 5.4 million phased 25. Cotterman,C.W. (1941) Relatives and human genetic analysis. Sci. human variants validated by genetic inheritance from sequencing a Mon., 53, 227–234. three-generation 17-member pedigree. Genome Res., 27, 157–164. 26. Li,X., Battle,A., Karczewski,K.J., Zappala,Z., Knowles,D.A., 38. Purcell,S., Neale,B., Todd-Brown,K., Thomas,L., Ferreira,M.A.R., Smith,K.S., Kukurba,K.R., Wu,E., Simon,N. and Montgomery,S.B. Bender,D., Maller,J., Sklar,P., de Bakker,P.I.W., Daly,M.J. et al. (2014) Transcriptome sequencing of a large human family identifies (2007) PLINK: a tool set for whole-genome association and the impact of rare noncoding variants. Am. J. Hum. Genet., 95, population-based linkage analyses. Am. J. Hum. Genet., 81, 559–575. 245–256. 39. R Core Team (2018) R: A Language and Environment for Statistical 27. Oda,H., Beck,D.B., Kuehn,H.S., Sampaio Moura,N., Hoffmann,P., Computing. Ibarra,M., Stoddard,J., Tsai,W.L., Gutierrez-Cruz,G., Gadina,M. 40. Galván-Femenı́a,I., Graffelman,J. and Barceló-I-Vidal,C. (2017) et al. (2019) Second case of HOIP deficiency expands clinical features Graphics for relatedness research. Mol. Ecol. Resour., 17, 1271–1282. and defines inflammatory transcriptome regulated by LUBAC. Front. 41. Staples,J., Qiao,D., Cho,M.H., Silverman,E.K., Nickerson,D.A., Immunol., 10, 479. Below,J.E. and University of Washington Center for Mendelian 28. Boutboul,D., Kuehn,H.S., Van de Wyngaert,Z., Niemela,J.E., Genomics (2014) PRIMUS: rapid reconstruction of pedigrees from Callebaut,I., Stoddard,J., Lenoir,C., Barlogis,V., Farnarier,C., Vely,F. genome-wide estimates of identity by descent. Am. J. Hum. Genet., et al. (2018) Dominant-negative IKZF1 mutations cause a T, B, and 95, 553–564. myeloid cell combined immunodeficiency. J. Clin. Invest., 128, 3071–3087.