Cell type-specific raw read counts and cis-eQTL summary statistics in healthy blood and gut

The CEDAR-2 dataset contains cell type-specific cis-eQTL summary statistics for all SNPs-gene pairs in: 27 distinct blood cell types sorted via FACS/MACS prior to bulk RNA sequencing (n = ~200 healthy individuals) 401 leaves/nodes in the hierchical tree constructed from single-cell RNA sequencing data of gut biopsies, including the terminal ileum, transverse colon, and rectum (n = 57 healthy individuals) This record contains: raw read counts after gene filtering in blood and gut permutational summary statistics of cis-eQTLs in blood and gut conditional summary statistics of cis-eQTLs in blood and gut nominal summary statistics of cis-eQTLs in blood The nominal summary statistics of cis-eQTLs in gut have been splitted in 14 other records. Raw read counts after gene filtering (blood): Demultiplexing and FASTQ conversion were performed using bcl2fastq (v2.20). Read quality was assessed with FastQC (v0.12.1) and multiQC (v0.9). Reads were mapped to the GRCh38 (Ensembl release 105) human genome build using STAR (v2.7.1a). The STAR re-implementation of the WASP algorithm was used with a custom VCF file containing both reference and alternative alleles for all SNPs. The alignments that did not pass WASP filtering or that overlapped indels were removed from the resulting BAM files using samtools (v1.9). Alignment metrics were collected using Picard CollectRnaSeqMetrics (v2.7.1) (STable 3). Matching of genomic and transcriptome genotypes was evaluated with QTLtools mbv (v1.3.1) (SFig. 2A). Unstranded gene counts were generated with HTSeq (v0.6.1p1). If samples were split across two sequencing lanes, we summed up the respective counts. To detect mislabeled samples, reads counts were normalized using the variance stabilizing transformation from the DESeq2 R package and a t-SNE analysis was conducted using the 500 most variable genes (SFig. 2). For each blood cell type, we filtered out the genes with less than 5 counts in more than 80% of samples. Raw read counts after gene filtering (gut): Raw sequencing data were preprocessed with the Cellranger software version 7.1.0 with standard parameters and GRCh38 (Ensembl release 103) human genome build as reference. Processed counts were further analyzed in R (version 4.3.1) within the Seurat tools ecosystem (Seurat version 4.1.3 [Hao et al., 2021]). For each sample, mRNA read counts and hashtag barcode read counts were loaded to R and quality-checked. We excluded (i) hashtag barcodes with < 300 reads per individual, (ii) cells with ≥ 50% mitochondrial reads, (iii) cells with < 200 different genes expressed, (iv) cells with < 5 hashtag barcode reads. Demultiplexing of cells was done using two different algorithms implemented in Seurat: HTODemux [Stoeckius et al., 2018] (positive.quantile 0.999, clusterization function - kmeans) and MULTIseqDemux [McGinnis et al., 2019] with automated threshold finding. Cells identified as singlets of the same tissue type by the two methods were kept for further analysis. To meet the RAM usage requirements, variable features selection was done in five sample batches with the “mean.var.plot" algorithm (3,000 features per batch). We retained the intersection between the five batches. We then integrated the sample batches in the space defined by the first 50 principal components using Harmony [Korsunsky et al., 2019]. The UMAP plots shown throughout the manuscript are in this coordinate system. Yet, to better differentiate cell types and improve the clustering, we further split the data in two stages. We first used Hashtag information to separate cells by anatomical location (IL, TC, RE). Within each anatomical location we repeated the variable feature selection and integration process. In each of these three sub-datasets, we identified cell clusters with the Louvain algorithm. Based on marker gene expression, we assigned clusters to three groups: immune (expressing PTPRC), epithelial (expressing EPCAM), and endothelial plus other cells (expressing VWF, PECAM1, CDH5 or not included in previous groups). Within each one of these nine sub-groups, we again repeated variable features selection, integration and clustering (Louvain clustering algorithm with resolution parameter 1.5). This yielded a total of 276 clusters across the nine data sets. We computed, for each of the 276 location- and cell-type specific cell clusters (obtained as described above), the mean coordinate vector in the space of the first 50 Harmony coordinates. Then we computed the Euclidean distance between these vectors followed by hierarchical clustering (function hclust, stats R package, “complete” agglomeration method). The final dendrogram was constructed with the dendextend [Galili, 2015] R package. We kept leaves and nodes in the hierarchical tree if the median cells per patient was ≥ 5 and that number of patients with cells in the leaf/node was ≥ 30. This left 401 leaves/nodes. Within each leaf or node, all cells were treated as pseudo-bulk, i.e., as if all reads were derived from one mega-cell. We kept genes if they were expressed in more than 10% of samples and with mean proportion of reads > 5e-7. Cis-eQTL analysis (blood): For each blood cell type, we normalized the raw read counts using the DESeq2 R package and residualized them for age, sex, RNA extraction method, proportion of reads in each sequencing batch, top 3 genotype principal components (PCs) and top 13-36 expression PCs to maximize the number of cis-eQTLs. STable 4 reports the number of samples, genes and top expression PCs for each cell type. We removed variants with MAF ≤ 0.05 in the retained samples using bcftools (v1.9) for each cell type. We performed eQTL mapping using QTLtools in a 2 Mb window centered at the transcription start site and with the integrated rank normal transformation of the phenotypes. The p-values were corrected for multiple testing within each window by permutation (10,000 permutations) and within each cell type using the false discovery rate (FDR). eQTL with “within cell type FDR” ≤ 0.05 were considered significant. Cis-eQTL analysis (gut): For each leaf and node in the hierarchical tree, gene expression data were normalized using DESeq2, residualized for age, sex and five genotype PCs, and corrected for hidden confounders utilizing the probabilistic estimation of expression residuals (PEER) algorithm [Stegle et al., 2010]. Association between gene expression and alternate allele dosage was conducted for each SNP in a 2Mb-window centered on the gene’s transcription start site using QTLtools [Delaneau et al., 2017]. Nominal p-values were corrected for the realization of multiple tests in this cis-window by permutation, yielding a window-adjusted p-value for one lead SNP for every gene x leaf/node combination. Window-adjusted lead SNP p-values for all genes in a leaf/node were jointly used to compute a q-value using Storey & Tibshirani [2003].