OSCR

Sibling chimerism among microglia in marmosets.

Code ↔ Paper

6 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 6 matches
  1. [1] § Methods › Donor-of-origin analysis and detection of host-sibling doublets (Dropulation) ↔ src/java/org/broadinstitute/dropseqrna/barnyard/digitalallelecounts/sampleassignment/multisample/DetectDoublets.java, lines 64–123 · score 0.69 · DetectDoublets, donor likelihood, AssignCellsToSamples, cell barcodes, Dropulation, threshold
  2. [2] § Methods › Latent factor analysis ↔ R/preprocessing.R, lines 3793–3878 · score 0.68 · Pearson residuals, variable features, glm, latent, SCT, regression
  3. [3] § Methods › Clustering of cells using independent component analysis ↔ R/integration.R, lines 3267–3410 · score 0.58 · nearest neighbor, expression matrix, algorithm, variable, libraries, cells
  4. [4] § Methods › Binomial generalized linear mixed-effects model analysis ↔ R/differential_expression.R, lines 455–538 · score 0.56 · binomial generalized linear, predictor, vector, identity, model, fraction
  5. [5] § Methods › Latent factor analysis ↔ R/vst.R, lines 5–109 · score 0.53 · Pearson residuals, glm, latent, regression, transform, variable
  6. [6] § Methods › Gene expression analysis of host and sibling meta cells ↔ src/python/src/dropseq/eqtl/normalize_tensorqtl_expression.py, lines 80–143 · score 0.52 · edgeR, fold change, sum, gene

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Java · 704 lines · 38 KB · MIT · 1 match

  1. /*
  2. * MIT License
  3. *
  4. * Copyright 2017 Broad Institute
  5. *
  6. * Permission is hereby granted, free of charge, to any person obtaining a copy
  7. * of this software and associated documentation files (the "Software"), to deal
  8. * in the Software without restriction, including without limitation the rights
  9. * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  10. * copies of the Software, and to permit persons to whom the Software is
  11. * furnished to do so, subject to the following conditions:
  12. *
  13. * The above copyright notice and this permission notice shall be included in all
  14. * copies or substantial portions of the Software.
  15. *
  16. * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  17. * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  18. * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  19. * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  20. * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  21. * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  22. * SOFTWARE.
  23. */
  24. package org.broadinstitute.dropseqrna.barnyard.digitalallelecounts.sampleassignment.multisample;
  26. import java.io.File;
  27. import java.io.PrintStream;
  28. import java.text.DecimalFormat;
  29. import java.util.ArrayList;
  30. import java.util.Arrays;
  31. import java.util.Collections;
  32. import java.util.Comparator;
  33. import java.util.HashMap;
  34. import java.util.HashSet;
  35. import java.util.List;
  36. import java.util.Map;
  37. import java.util.Set;
  39. import org.apache.commons.lang3.StringUtils;
  40. import org.broadinstitute.barclay.argparser.Argument;
  41. import org.broadinstitute.barclay.argparser.CommandLineProgramProperties;
  42. import org.broadinstitute.dropseqrna.barnyard.GeneFunctionCommandLineBase;
  43. import org.broadinstitute.dropseqrna.barnyard.ParseBarcodeFile;
  44. import org.broadinstitute.dropseqrna.barnyard.digitalallelecounts.SNPUMIBasePileupIterator;
  45. import org.broadinstitute.dropseqrna.barnyard.digitalallelecounts.SortOrder;
  46. import org.broadinstitute.dropseqrna.barnyard.digitalallelecounts.sampleassignment.CellAssignmentUtils;
  47. import org.broadinstitute.dropseqrna.barnyard.digitalallelecounts.sampleassignment.CellCollectionSampleLikelihoodCollection;
  48. import org.broadinstitute.dropseqrna.barnyard.digitalallelecounts.sampleassignment.SampleGenotypeProbabilities;
  49. import org.broadinstitute.dropseqrna.barnyard.digitalallelecounts.sampleassignment.SampleGenotypeProbabilitiesIterator;
  50. import org.broadinstitute.dropseqrna.cmdline.CustomCommandLineValidationHelper;
  51. import org.broadinstitute.dropseqrna.cmdline.DropSeq;
  52. import org.broadinstitute.dropseqrna.utils.AssertSequenceDictionaryIntersection;
  53. import org.broadinstitute.dropseqrna.utils.FileListParsingUtils;
  54. import org.broadinstitute.dropseqrna.utils.FileUtils;
  55. import org.broadinstitute.dropseqrna.utils.GroupingIterator;
  56. import org.broadinstitute.dropseqrna.utils.VCFUtils;
  57. import org.broadinstitute.dropseqrna.utils.io.ErrorCheckingPrintStream;
  58. import org.broadinstitute.dropseqrna.utils.readiterators.IgnoreGeneAnnotationTagger;
  59. import org.broadinstitute.dropseqrna.utils.readiterators.PCRDuplicateFilteringIterator;
  60. import org.broadinstitute.dropseqrna.utils.readiterators.SamFileMergeUtil;
  61. import org.broadinstitute.dropseqrna.utils.readiterators.SamHeaderAndIterator;
  62. import org.broadinstitute.dropseqrna.vcftools.SampleAssignmentVCFUtils;
  64. import htsjdk.samtools.SAMRecord;
  65. import htsjdk.samtools.SAMSequenceDictionary;
  66. import htsjdk.samtools.SamReaderFactory;
  67. import htsjdk.samtools.util.CloseableIterator;
  68. import htsjdk.samtools.util.IOUtil;
  69. import htsjdk.samtools.util.Interval;
  70. import htsjdk.samtools.util.IntervalList;
  71. import htsjdk.samtools.util.Log;
  72. import htsjdk.samtools.util.PeekableIterator;
  73. import htsjdk.variant.variantcontext.VariantContext;
  74. import htsjdk.variant.vcf.VCFFileReader;
  75. import picard.cmdline.StandardOptionDefinitions;
  76. import picard.nio.PicardHtsPath;
  78. @CommandLineProgramProperties(summary = "Detect Doublets in Dropulation Data. Uses the outputs of AssignCellsToSamples to make decisions. It's highly recommended to use the VCF output from AssignCellsToSamples as input, as"
  79. + "the memory usage of a full VCF may be prohibitive compared to the AssignCellsToSamples VCF, which contains only the variants that were observed in the data. This also greatly speeds up"
  80. + "analysis.", oneLineSummary = "Detect Doublets in Dropulation Data", programGroup = DropSeq.class)
  81. public class DetectDoublets extends GeneFunctionCommandLineBase {
  83. private static final Log log = Log.getInstance(DetectDoublets.class);
  85. @Argument(shortName = StandardOptionDefinitions.INPUT_SHORT_NAME, doc = "The input SAM or BAM file to analyze. This argument can accept wildcards, or a file with the suffix .bam_list that contains the locations of multiple BAM files", minElements = 1)
  86. public List<PicardHtsPath> INPUT_BAM;
  88. @Argument(doc = "The input VCF file to analyze. Use the output VCF from AssignCellsToSamples to save memory.")
  89. public PicardHtsPath VCF;
  91. @Argument(doc = "The output likelihood file from AssignCellsToSamples")
  92. public File SINGLE_DONOR_LIKELIHOOD_FILE;
  94. @Argument(shortName = StandardOptionDefinitions.OUTPUT_SHORT_NAME, doc = "Output file of doublet likelihoods. This supports zipped formats like gz and bz2.")
  95. public File OUTPUT;
  97. @Argument(doc = "Output file of per-sample-pair doublet likelihoods. Insted of just the best pair as seen in the OUTPUT file, "
  98. + "this outputs all tested pairs for each cell. This supports zipped formats like gz and bz2.", optional = true)
  99. public File OUTPUT_ALL_PAIRS = null;
  101. @Argument(doc = "Output file of per-snp/sample-pair doublet likelihoods. Insted of just the best pair as seen in the OUTPUT file, "
  102. + "this outputs all tested pairs for each cell, for each SNP. This file can get pretty obnoxiously huge."
  103. + "This supports zipped formats like gz and bz2.", optional = true)
  104. public File OUTPUT_PER_SNP = null;
  106. @Argument(doc = "The cell barcode tag. If there are no reads with this tag, the program will assume that all reads belong to the same cell and process in single sample mode.")
  107. public String CELL_BARCODE_TAG = "XC";
  109. @Argument(doc = "The molecular barcode tag.")
  110. public String MOLECULAR_BARCODE_TAG = "XM";
  112. @Argument(doc = "The edit distance that molecular barcodes should be combined at within a gene/SNP.")
  113. public Integer EDIT_DISTANCE = 1;
  115. @Argument(doc = "The map quality of the read to be included.")
  116. public Integer READ_MQ = 10;
  118. @Argument(doc = "Override NUM_CORE_BARCODES and process reads that have the cell barcodes in this file instead. The file has 1 column with no header.", optional = false)
  119. public File CELL_BC_FILE = null;
  121. @Argument(doc = "The minimum genotype quality for a variant. Set this value to 0 to not filter by GQ scores if they are present, or to -1 to completely "
  122. + "ignore GQ values if they are not set in the genotype info field. If the GQ field is not set in the VCF header, this will be set to -1 by default.")
  123. public Integer GQ_THRESHOLD = 30;
  125. @Argument(doc = "A file with a list of samples in the VCF to consider as samples in the doublets. This subsets the VCF into a smaller data set containing only the samples listed. "
  126. + "The file has 1 column with no header. If this list contains only one donor and no contaminating donors were found via single donor assignment, "
  127. + "doublet detection calculations will not take place. Instead, the program will emit a default output for each cell, and the program "
  128. + "will then quit with an exit status.", optional = false)
  129. public File SAMPLE_FILE;
  131. @Argument(doc = "Instead of useing base qualities to determine error rate, use a fixed error rate instead. This is rounded to the nearest phread score internally.", optional = true)
  132. public Double FIXED_ERROR_RATE = null;
  134. @Argument(doc = "Caps the base error rate at a maximum probability so no SNP can be weighed more than this value. For example, if this value was 0.01, "
  135. + "then a base quality 30 value (normally an erro rate of 0.001) would become 0.01. With the same threshold, a base with an error rate of 0.1 would be unaffected.", optional = true)
  136. public Double MAX_ERROR_RATE = null;
  138. @Argument(doc = "A file that contains an estimate of how much ambient RNA is in each cell. This is a fractional estimate between 0 and 1. File is tab seperated, with 2 columns:"
  139. + "cell_barcode and frac_contamination. When supplied along side the ALLELE_FREQUENCY_ESTIMATE_FILE, this modifies the likelihood error rates to take into account how often"
  140. + "the allele observed can be drawn from ambient RNA. We use cellbender remove background [https://github.com/broadinstitute/CellBender] to estimate the "
  141. + "number of transcripts before and after ambient cleanup to define the fraction of transcripts that come from ambient RNA.", optional = true)
  142. public File CELL_CONTAMINATION_ESTIMATE_FILE = null;
  144. @Argument(doc = "A file that contains an estimate of the allele frequency expected for each SNP across donors. The best estimate of this will come from the fraction of reference and alternate allele"
  145. + "UMIs that are observed at each snp site. This report can be generated via GatherDigitalAlleleCounts. This is a fractional estimate between 0 and 1. File is tab seperated, with at least 3 columns:"
  146. + "chromosome, position, maf_umi. When supplied and CELL_CONTAMINATION_ESTIMATE_FILE is provided, this modifies the likelihood error rates to take into account how often "
  147. + "the allele observed can be drawn from ambient RNA.", optional = true)
  148. public File ALLELE_FREQUENCY_ESTIMATE_FILE = null;
  150. @Argument(doc = "At least <FRACTION_SAMPLES_PASSING> samples must have genotype scores >= GQ_THRESHOLD for the variant in the VCF to be included in the analysis.")
  151. public double FRACTION_SAMPLES_PASSING = 0.5;
  153. @Argument(doc = "A list of chromosomes to omit from the analysis. The default is to omit the sex chromosomes.")
  154. public List<String> IGNORED_CHROMOSOMES = new ArrayList<>(Arrays.asList("X", "Y", "MT"));
  156. /**
  157. * this produces worse results with in-silico testing, generating more misclassifications of singlets as doublets.
  158. *
  159. * @Argument(doc="For SNPs that don't have a high quality genotype in the VCF, should we infer a global penalty per SNP
  160. * and apply it to donors for that SNP that are not confidently called?") public boolean
  161. * USE_MISSING_DATA=true;
  162. */
  163. private boolean USE_MISSING_DATA = false;
  165. @Argument(doc = "Force evaluation of the doublet at the given mixture ratio. Should be a number between 0 and 1.", optional = true)
  166. public Double FORCED_RATIO = 0.8;
  168. @Argument(doc = "Should cells that were assigned to a donor not on the sample list be tested? When enabled, "
  169. + "this forces the program to load genotype information for all donors seen at least once, even when not on the sample list. If there are many off-target assignments,"
  170. + "this can use large amounts of memory. Set to false to skip testing cells that you'll probably discard as incorrectly assigned later anyway, set to true to "
  171. + "test each cell not on the donor lists against all possible donors on the list.")
  172. public Boolean TEST_UNEXPECTED_DONORS = true;
  174. @Argument(doc = "For each cell, when comparing donor pairs to each other, scale the likelihoods to the number of UMIs. "
  175. + "This provides an additional penalty score to donor pairs with more incomplete information that makes donor pairs more comparable.")
  176. public Boolean SCALE_LIKELIHOODS = true;
  178. @Argument(doc = "EXPERIMENTAL!!! Run the program in DNA Mode. In this mode, reads should have a cell barcode, but will be missing gene annotations and UMIs. All reads will be "
  179. + "accepted as passing, and each read (or read pair) will be treated as a single UMI If the data is PCR Duplicate marked, duplicate reads will be filtered. ")
  180. public Boolean DNA_MODE = false;
  182. @Argument(doc = "Set to false exclude donor cells that are unable to be tested for doublets.")
  183. public Boolean WRITE_ALL_DONOR_BARCODES = true;
  185. @Argument(doc = "The value to write for the best pair p-value for untested cells.")
  186. public Double MISSING_BEST_PAIR_PVALUE = 1.0E-101d;
  188. /*
  189. * @Argument
  190. * (doc="The model tests the best donor against the all the possible second most likely donors to find the pair that best explain the data. Sort the donors by their single donor likelihood score,"
  191. * +
  192. * "and only test the best <MIN_DONOR_PAIRS_TESTED> donors. This can reduce the number of tests / runtime / memory for large pools while producing approximately the same output."
  193. * )
  194. */
  196. private Integer MIN_DONOR_PAIRS_TESTED = Integer.MAX_VALUE;
  198. // @Argument(doc="See MIN_DONOR_PAIRS_TESTED. This limits the number of donors tested to a fraction of the total number
  199. // of donors. The number of donors used is the maximum of FRACTION_DONOR_PAIRS_TESTED and MIN_DONOR_PAIRS_TESTED.")
  201. private Double FRACTION_DONOR_PAIRS_TESTED = 1.0;
  203. private final String SNP_TAG = "YS";
  205. private static DecimalFormat mixtureFormat = new DecimalFormat("#.###");
  207. @Override
  208. protected int doWork() {
  209. if (CELL_CONTAMINATION_ESTIMATE_FILE != null) {
  210. IOUtil.assertFileIsReadable(CELL_CONTAMINATION_ESTIMATE_FILE);
  211. }
  212. if (ALLELE_FREQUENCY_ESTIMATE_FILE != null) {
  213. IOUtil.assertFileIsReadable(ALLELE_FREQUENCY_ESTIMATE_FILE);
  214. }
  216. List<String> donorList = ParseBarcodeFile.readCellBarcodeFile(this.SAMPLE_FILE);
  217. log.info("Number of donors in donor list [" + donorList.size() + "]");
  218. int pairsToTest = getNumberOfDonorsToTest(donorList.size(), this.MIN_DONOR_PAIRS_TESTED, this.FRACTION_DONOR_PAIRS_TESTED);
  220. PrintStream perDonorWriter = null;
  221. if (OUTPUT_ALL_PAIRS != null) {
  222. perDonorWriter = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(this.OUTPUT_ALL_PAIRS));
  223. writeAssignmentHeader(perDonorWriter, pairsToTest, false, true);
  224. }
  226. PrintStream perSNPWriter = null;
  227. if (OUTPUT_PER_SNP != null) {
  228. IOUtil.assertFileIsWritable(this.OUTPUT_PER_SNP);
  229. perSNPWriter = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(this.OUTPUT_PER_SNP));
  230. writePerSNPHeader(perSNPWriter);
  231. }
  233. PrintStream writer = new ErrorCheckingPrintStream(IOUtil.openFileForWriting(this.OUTPUT));
  234. writeAssignmentHeader(writer, pairsToTest, true, false);
  236. final VCFFileReader vcfReader = new VCFFileReader(this.VCF.toPath(), false);
  237. final SamHeaderAndIterator headerAndIter = SamFileMergeUtil.mergeInputPaths(
  238. PicardHtsPath.toPaths(this.INPUT_BAM), false, SamReaderFactory.makeDefault());
  239. AssertSequenceDictionaryIntersection.assertIntersectionObjectVcf(
  240. headerAndIter.header, "BAM INPUT(S)", this.VCF.toPath(), log);
  242. // extract the sequence dictionary to build the interval list.
  243. SAMSequenceDictionary vcfDict = vcfReader.getFileHeader().getSequenceDictionary();
  245. // disable GQ filter if it's not in the header.
  246. if (!VCFUtils.GQInHeader(vcfReader)) {
  247. this.GQ_THRESHOLD = -1;
  248. log.info("Genotype Quality [GQ] not found in header. Disabling GQ_THRESHOLD parameter");
  249. }
  251. CellCollectionSampleLikelihoodCollection cslc = CellCollectionSampleLikelihoodCollection.parseFile(this.SINGLE_DONOR_LIKELIHOOD_FILE);
  253. // A map where the key is the cell barcode, and the value is the best donor.
  254. Map<String, String> bestDonorForCell = getBestDonorForCell(this.SINGLE_DONOR_LIKELIHOOD_FILE);
  256. // Filter the single donor assignment by the donor list if requested.
  257. if (!this.TEST_UNEXPECTED_DONORS)
  258. bestDonorForCell = filterDonorMap(bestDonorForCell, donorList);
  260. // read in the per-cell penalty score and validate against the best donor per cell to make sure you have penalties for
  261. // every donor.
  262. Map<String, Double> contaminationMap = CellAssignmentUtils.getCellContamination(this.CELL_CONTAMINATION_ESTIMATE_FILE, bestDonorForCell.keySet());
  263. Map<Interval, Double> variantMinorAlleleFrequency = CellAssignmentUtils.getMinorAlleleFrequencyMap(this.ALLELE_FREQUENCY_ESTIMATE_FILE);
  265. // all donors is the input set of donors plus any best donor for a cell.
  266. // this list of donors is restricted to the sample list if TEST_UNEXPECTED_DONORS=false.
  267. Set<String> allDonors = new HashSet<>();
  268. allDonors.addAll(donorList);
  269. allDonors.addAll(new HashSet<>(bestDonorForCell.values()));
  270. List<String> allDonorsList = new ArrayList<>(allDonors);
  272. log.info("Number of donors that can be either donor in a donor pair [" + allDonorsList.size() + "]");
  274. // Keep all the barcodes that are in the barcode list AND in the single donor assignments.
  275. List<String> cellBarcodes = getCellBarcodes(this.CELL_BC_FILE, bestDonorForCell, this.TEST_UNEXPECTED_DONORS);
  277. // If there is only one donor at this point, doublet detection should not continue.
  278. if (allDonorsList.size()<2) {
  279. singleDonorGracefulExit(bestDonorForCell, writer, perDonorWriter, perSNPWriter);
  280. return 0;
  281. }
  283. // Pass a list of all donors requested + best calls if TEST_UNEXPECTED_DONORS is true.
  284. // set the %passing to be 0, since the snpIntervals will properly filter on the right set of donors, but you want a
  285. // super-set of donors available
  286. // so if you call donors A-D, but the single cell assigned E, then E would still be an available donor in the genotype
  287. // matrix.
  288. PeekableIterator<VariantContext> vcfIterator = SampleAssignmentVCFUtils.getVCFIterator(vcfReader, allDonorsList, false, this.GQ_THRESHOLD,
  289. this.FRACTION_SAMPLES_PASSING, IGNORED_CHROMOSOMES, log);
  291. GenotypeMatrix genotypeMatrix = new GenotypeMatrix(vcfIterator, this.GQ_THRESHOLD, allDonorsList);
  292. vcfIterator.close();
  294. Map<Interval, Double> genotypeQuality = genotypeMatrix.getAverageGenotypeQuality();
  296. final IntervalList snpIntervals = new IntervalList(vcfDict);
  297. snpIntervals.addall(genotypeMatrix.getSNPIntervals());
  298. // a requested early exit if there are no SNPs.
  299. if (snpIntervals.getIntervals().isEmpty()) {
  300. log.error("No SNP intervals detected! Check to see if your VCF filter thresholds are too restrictive!");
  301. return 1;
  302. }
  304. PeekableIterator<List<SampleGenotypeProbabilities>> sampleGenotypeIterator = prepareIterator(snpIntervals, cellBarcodes, genotypeQuality);
  306. int cellCount = 0;
  307. int reportInterval = 100;
  308. final Set<String> writtenBarcodes = new HashSet<>();
  309. log.info("Calling doublets");
  310. if (!sampleGenotypeIterator.hasNext()) {
  311. log.warn("No Cells found for analysis.");
  312. } else {
  314. while (sampleGenotypeIterator.hasNext()) {
  315. cellCount++;
  316. if (cellCount % reportInterval == 0)
  317. log.info("Tested cell #" + cellCount);
  318. List<SampleGenotypeProbabilities> probs = sampleGenotypeIterator.next();
  320. String cell = probs.get(0).getCell();
  321. String bestDonor = bestDonorForCell.get(cell);
  322. if (bestDonor == null)
  323. throw new IllegalStateException("Cell [" + cell + "] has no best donor assignment in file.");
  325. VariantDataFactory f = null;
  326. f = new VariantDataFactory(cell, probs, genotypeMatrix, FIXED_ERROR_RATE, USE_MISSING_DATA, MAX_ERROR_RATE, contaminationMap,
  327. variantMinorAlleleFrequency);
  329. FindOptimalDonorMixture fodm = new FindOptimalDonorMixture(f);
  331. // AllPairedSampleAssignmentsForCell allAssignments = fodm.findBestDonorPair(bestDonor, donorList, FORCED_RATIO);
  332. // List<String> donorsThisCell = getExpectedSecondDonorsRankedByLikelihood(cell, cslc, pairsToTest, donorList, bestDonor);
  333. AllPairedSampleAssignmentsForCell allAssignments = fodm.findBestDonorPair(bestDonor, donorList, FORCED_RATIO, SCALE_LIKELIHOODS);
  334. SamplePairAssignmentForCell best = allAssignments.getBestAssignment();
  336. // edge case: assignment is null because there's no data for this cell. This only happens
  337. // when a user has a cell selection error or similar and attempts to call cell barcodes that aren't cells.
  338. if (best==null) {
  339. log.warn("No best pair found for cell ["+cell+"] due to no informative UMIs. Cell selection or similar problem?");
  340. best= SamplePairAssignmentForCell.constructEmptyResult(cell, bestDonorForCell.get(cell));
  341. }
  343. double bestPairPvalue = allAssignments.getBestPairPvalue();
  344. writeAssignment(best, bestPairPvalue, writer, false);
  345. if (OUTPUT_ALL_PAIRS != null) {
  346. writeAssignment(allAssignments.getBestAssignment(), null, perDonorWriter, true);
  347. // apply ordering to other assignments for stability of outputs. Sort by 2nd donor name.
  348. Comparator<SamplePairAssignmentForCell> comparator = java.util.Comparator.comparing(SamplePairAssignmentForCell::getSampleTwo,
  349. java.util.Comparator.naturalOrder());
  350. List<SamplePairAssignmentForCell> allOther = allAssignments.getOtherAssignments();
  351. Collections.sort(allOther, comparator);
  352. for (SamplePairAssignmentForCell other : allOther)
  353. writeAssignment(other, null, perDonorWriter, true);
  354. }
  355. reportResultsPerSNP(cell, f, bestDonor, donorList, allAssignments, perSNPWriter);
  356. writtenBarcodes.add(cell);
  357. }
  358. }
  360. if (WRITE_ALL_DONOR_BARCODES) {
  361. final List<String> remainingCellBarcodes = new ArrayList<>(bestDonorForCell.keySet());
  362. remainingCellBarcodes.removeAll(writtenBarcodes);
  363. Collections.sort(remainingCellBarcodes);
  364. for (final String cell : remainingCellBarcodes) {
  365. log.warn("No best pair found for cell [" + cell + "] due to no informative UMIs. Cell selection or similar problem?");
  366. final SamplePairAssignmentForCell best = SamplePairAssignmentForCell.constructEmptyResult(cell, bestDonorForCell.get(cell));
  367. writeAssignment(best, MISSING_BEST_PAIR_PVALUE, writer, false);
  368. }
  369. }
  371. if (OUTPUT_PER_SNP != null)
  372. perSNPWriter.close();
  373. if (OUTPUT_ALL_PAIRS != null)
  374. perDonorWriter.close();
  375. writer.close();
  377. log.info("Finished!");
  378. return 0;
  379. }
  382. /**
  383. * In the strange edge case where there is only a single donor to be tested, write out a default output file instead of going through testing.
  384. * This function additionally closes all potentially open writers and runs logging.
  385. * @param bestDonorForCell A map containing cell barcodes and the best donor for each cell.
  386. * @param writer The file to write to per-cell outputs.
  387. * @param perDonorWriter Closes this file if not null.
  388. * @param perSNPWriter Closes this file if not null.
  389. */
  390. void singleDonorGracefulExit(Map<String, String> bestDonorForCell, PrintStream writer,
  391. PrintStream perDonorWriter, PrintStream perSNPWriter) {
  392. // clean up more detailed file writers.
  393. if (OUTPUT_ALL_PAIRS!=null) perDonorWriter.close();
  394. if (OUTPUT_PER_SNP != null) perSNPWriter.close();
  395. // write a default output per cell close results and quit.
  396. log.error("The donor file only contained a single donor, and no additional donors were detected by single donor assignment. Doublet detection will not continue. "
  397. + "A default output will be written to perserve downstream pipeline functionality.");
  398. writeSingleDonorEdgeCaseOutput(bestDonorForCell, writer);
  399. }
  401. /**
  402. * In the strange edge case where there is only a single donor to be tested, write out a default output file instead of going through testing.
  403. * @param bestDonorForCell A map containing cell barcodes and the best donor for each cell.
  404. * @param writer The file to write to
  405. */
  406. void writeSingleDonorEdgeCaseOutput(Map<String, String> bestDonorForCell, PrintStream writer) {
  407. for (String cell: bestDonorForCell.keySet()) {
  408. SamplePairAssignmentForCell best = SamplePairAssignmentForCell.constructEmptyResult(cell, bestDonorForCell.get(cell));
  409. writeAssignment(best, MISSING_BEST_PAIR_PVALUE, writer, false);
  410. }
  411. writer.close();
  413. }
  415. /**
  416. * Filter the map of cell barcode -> donor to only retain cell barcodes where the assigned donor is in the donor list.
  417. *
  418. * @param bestDonorForCell
  419. * @param donorList
  420. * @return A submap of the input map where all values are contained in the donor list.
  421. */
  422. Map<String, String> filterDonorMap(Map<String, String> bestDonorForCell, List<String> donorList) {
  423. Set<String> dl = new HashSet<String>(donorList);
  424. Map<String, String> result = new HashMap<String, String>();
  426. for (String cellBC : bestDonorForCell.keySet()) {
  427. String donor = bestDonorForCell.get(cellBC);
  428. if (dl.contains(donor))
  429. result.put(cellBC, donor);
  430. }
  431. return result;
  432. }
  434. private void reportResultsPerSNP(final String cell, final VariantDataFactory variantFactory, final String bestDonor, final List<String> vcfSamples,
  435. final AllPairedSampleAssignmentsForCell allAssignments, final PrintStream out) {
  436. if (out == null)
  437. return;
  438. List<String> other = FindOptimalDonorMixture.getNonPrimarySamples(bestDonor, vcfSamples);
  439. for (String o : other) {
  440. VariantDataCollection vdc = variantFactory.getVariantData(bestDonor, o);
  441. SamplePairAssignmentForCell mixtureResult = allAssignments.getAssignmentForDonorPair(bestDonor, o);
  442. List<VariantData> vdList = vdc.getVariantData();
  443. double mixture = mixtureResult.getMixture();
  444. for (VariantData vd : vdList)
  445. writePerSNPReport(cell, vd, bestDonor, o, mixture, out);
  446. }
  447. }
  449. private void writePerSNPReport(final String cell, final VariantData vd, final String sampleOne, final String sampleTwo, final Double mixture,
  450. final PrintStream out) {
  451. /*
  452. * if (mixture==null) { String [] line = {cell, sampleOne, sampleTwo, vd.getSNPInterval().getContig(),
  453. * Integer.toString(vd.getSNPInterval().getStart()), vd.getGenotypeOne().toString(), vd.getGenotypeTwo().toString(),
  454. * Integer.toString(vd.getGenotypeCountReference()), Integer.toString(vd.getGenotypeCountAlternate()), "NA",
  455. * Double.toString(vd.getLogLikelihood(1)), Double.toString(vd.getLogLikelihood(0))}; String h = StringUtils.join(line,
  456. * "\t"); out.println(h); return; }
  457. */
  459. String[] line = { cell, sampleOne, sampleTwo, vd.getSNPInterval().getContig(), Integer.toString(vd.getSNPInterval().getStart()),
  460. vd.getGenotypeOne().toString(), vd.getGenotypeTwo().toString(), Integer.toString(vd.getGenotypeCountReference()),
  461. Integer.toString(vd.getGenotypeCountAlternate()), Double.toString(vd.getLogLikelihood(mixture)), Double.toString(vd.getLogLikelihood(1)),
  462. Double.toString(vd.getLogLikelihood(0)) };
  463. String h = StringUtils.join(line, "\t");
  464. out.println(h);
  465. }
  467. private void writePerSNPHeader(final PrintStream out) {
  468. String[] line = { "cell", "sampleOne", "sampleTwo", "chr", "pos", "genotype_S1", "genotype_S2", "refAlleleCount", "altAlleleCount",
  469. "likelihood_mixture", "likelihood_S1", "likelihood_S2" };
  470. String h = StringUtils.join(line, "\t");
  471. out.println(h);
  472. }
  474. private String convertNullToString(final Double x) {
  475. if (x == null)
  476. return ("NA");
  477. return Double.toString(x);
  478. }
  480. private void writeAssignmentHeader(final PrintStream out, final int pairsToTest, final boolean outputBestPairPvalue, final boolean writeScaledLikelihoods) {
  481. final List<String> paths = FileUtils.toAbsoluteStrings(PicardHtsPath.toPaths(this.INPUT_BAM));
  482. String bamList = StringUtils.join(paths, ",");
  484. final List<String> header = new ArrayList<>(Arrays.asList(
  485. "#INPUT_BAM=" + bamList, "INPUT_VCF=" + FileUtils.toAbsoluteString(this.VCF.toPath()),
  486. "DONOR_FILE=" + this.SAMPLE_FILE, "CELL_BC_FILE=" + CELL_BC_FILE, "GQ_THRESHOLD=" + Integer.toString(this.GQ_THRESHOLD),
  487. "FRACTION_SAMPLES_PASSING=" + Double.toString(this.FRACTION_SAMPLES_PASSING), "FORCED_RATIO=" + convertNullToString(FORCED_RATIO),
  488. "USE_MISSING_DATA=" + USE_MISSING_DATA, "READ_MQ=" + Integer.toString(this.READ_MQ),
  489. "FIXED_ERROR_RATE=" + convertNullToString(this.FIXED_ERROR_RATE), "MAX_ERROR_RATE=" + convertNullToString(this.MAX_ERROR_RATE),
  490. "LOCUS_FUNCTION=" + this.LOCUS_FUNCTION_LIST.toString(), "PAIRS_TO_TEST=" + Integer.toString(pairsToTest)));
  492. if (this.CELL_CONTAMINATION_ESTIMATE_FILE != null && this.ALLELE_FREQUENCY_ESTIMATE_FILE != null) {
  493. header.add("CELL_CONTAMINATION_ESTIMATE_FILE=" + CELL_CONTAMINATION_ESTIMATE_FILE.getAbsolutePath());
  494. header.add("ALLELE_FREQUENCY_ESTIMATE_FILE=" + ALLELE_FREQUENCY_ESTIMATE_FILE.getAbsolutePath());
  495. }
  497. String h = StringUtils.join(header, "\t");
  498. out.println(h);
  500. writeAssignmentColumnNames(out, outputBestPairPvalue, writeScaledLikelihoods);
  501. }
  503. public static void writeAssignmentColumnNames(final PrintStream out, final boolean outputBestPairPvalue, final boolean writeScaledLikelihoods) {
  504. List<String> line = new ArrayList<String>(Arrays.asList("cell", "sampleOneMixtureRatio", "sampleOne", "sampleOneLikelihood", "sampleTwo",
  505. "sampleTwoLikelihood", "mixedSample", "mixedSampleLikelihood", "num_paired_snps", "num_inform_snps", "num_umi", "num_inform_umis",
  506. "lr_test_stat", "sampleOneWrongAlleleCount", "num_homozygous_inform_umis_s1",
  507. "sampleTwoWrongAlleleCount", "num_homozygous_inform_umis_s2", "bestLikelihood", "bestSample", "doublet_pval"));
  509. if (outputBestPairPvalue)
  510. line.add("best_pair_pvalue");
  512. if (writeScaledLikelihoods) {
  513. line.add("bestLikelihoodScaled");
  514. }
  516. String header = StringUtils.join(line, "\t");
  517. out.println(header);
  518. }
  520. public static void writeAssignment(final SamplePairAssignmentForCell assignment, final Double bestPairPvalue, final PrintStream out,
  521. final boolean writeScaledLikelihood) {
  522. String mixture = mixtureFormat.format(assignment.getMixture());
  524. List<String> line = new ArrayList<String>(
  525. Arrays.asList(assignment.getCellBarcode(), mixture, assignment.getSampleOne(), Double.toString(assignment.getSampleOneSingleLikelihood()),
  526. assignment.getSampleTwo(), Double.toString(assignment.getSampleTwoSingleLikelihood()), assignment.getCombinedDonorName(),
  527. Double.toString(assignment.getDoubletLikelihood()), Integer.toString(assignment.getNumSNPs()),
  528. Integer.toString(assignment.getNumInformativeSNPs()), Integer.toString(assignment.getNumUMIs()),
  529. Integer.toString(assignment.getNumInformativeUMIs()), Double.toString(assignment.getDoubletLikelihoodRatio()),
  530. Integer.toString(assignment.getImpossibleAllelesSampleOne()), Integer.toString(assignment.getNumInformativeHomozygousUMIsSampleOne()),
  531. Integer.toString(assignment.getImpossibleAllelesSampleTwo()), Integer.toString(assignment.getNumInformativeHomozygousUMIsSampleTwo()),
  532. Double.toString(assignment.getBestLikelihood()), assignment.getBestSample(), Double.toString(assignment.getDoubletPvalue())));
  534. if (bestPairPvalue != null)
  535. line.add(bestPairPvalue.toString());
  537. if (writeScaledLikelihood)
  538. line.add(Double.toString(assignment.getScaledBestLikelihood()));
  540. String h = StringUtils.join(line, "\t");
  541. out.println(h);
  542. }
  544. public PeekableIterator<List<SampleGenotypeProbabilities>> prepareIterator(final IntervalList snpIntervals, List<String> cellBarcodes, Map<Interval, Double> genotypeQuality) {
  546. SamReaderFactory factory = SamReaderFactory.makeDefault().enable(SamReaderFactory.Option.EAGERLY_DECODE);
  547. SamHeaderAndIterator headerAndIter =
  548. SamFileMergeUtil.mergeInputPaths(PicardHtsPath.toPaths(this.INPUT_BAM), false, factory);
  550. // override the normal gene annotations with new ones before any other operations.
  551. // filter out PCR duplicates.
  552. if (this.DNA_MODE) {
  553. // Replace UMI tags with read names, set strand tag to read tag, set gene function to an accepted on. Overwrites ALL
  554. // reads tags.
  555. IgnoreGeneAnnotationTagger tagger = new IgnoreGeneAnnotationTagger(headerAndIter.iterator, this.GENE_NAME_TAG, this.GENE_STRAND_TAG,
  556. this.GENE_FUNCTION_TAG, this.LOCUS_FUNCTION_LIST, false, this.MOLECULAR_BARCODE_TAG, true);
  557. headerAndIter = new SamHeaderAndIterator(headerAndIter.header, (CloseableIterator<SAMRecord>) tagger.iterator());
  558. final PCRDuplicateFilteringIterator pcrDuplicateFilteringIterator = new PCRDuplicateFilteringIterator(headerAndIter.iterator);
  559. headerAndIter = new SamHeaderAndIterator(headerAndIter.header, pcrDuplicateFilteringIterator);
  560. }
  562. SNPUMIBasePileupIterator sbpi = new SNPUMIBasePileupIterator(headerAndIter, snpIntervals, GENE_NAME_TAG, GENE_STRAND_TAG, GENE_FUNCTION_TAG,
  563. LOCUS_FUNCTION_LIST, STRAND_STRATEGY, this.FUNCTIONAL_STRATEGY, this.CELL_BARCODE_TAG, this.MOLECULAR_BARCODE_TAG, this.SNP_TAG,
  564. GeneFunctionCommandLineBase.DEFAULT_FUNCTION_TAG, this.READ_MQ, false, cellBarcodes, genotypeQuality, SortOrder.CELL_SNP);
  566. final SAMSequenceDictionary dict = snpIntervals.getHeader().getSequenceDictionary();
  568. // gets a SampleGenotypeProbabilities for each cell.
  569. SampleGenotypeProbabilitiesIterator result = new SampleGenotypeProbabilitiesIterator(sbpi, dict, this.EDIT_DISTANCE, SortOrder.CELL_SNP);
  571. // clusters SampleGenotypeProbabilities objects across all cells for a SNP
  572. GroupingIterator<SampleGenotypeProbabilities> groupingIterator = new GroupingIterator<>(result, new Comparator<SampleGenotypeProbabilities>() {
  573. @Override
  574. public int compare(final SampleGenotypeProbabilities o1, final SampleGenotypeProbabilities o2) {
  575. int cmp = o1.getCell().compareTo(o2.getCell());
  576. return cmp;
  577. }
  578. });
  580. PeekableIterator<List<SampleGenotypeProbabilities>> peekableIter = new PeekableIterator<>(groupingIterator);
  581. return (peekableIter);
  582. }
  584. List<String> getCellBarcodes(File cellBarcodeFile, Map<String, String> bestDonorForCell, boolean testUnexpectedDonors) {
  585. List<String> cellBarcodes = ParseBarcodeFile.readCellBarcodeFile(this.CELL_BC_FILE);
  586. int numBarcodes = cellBarcodes.size();
  588. log.info("Number of cell barcodes in input file [" + numBarcodes + "]");
  590. if (!testUnexpectedDonors) {
  591. cellBarcodes.retainAll(bestDonorForCell.keySet());
  592. log.info("Number of cell barcodes after filtering to expected donors [" + cellBarcodes.size() + "]");
  593. double fracRemoved = (double) (numBarcodes - cellBarcodes.size()) / (double) numBarcodes;
  594. log.info("% cell barcodes not assigned to an expected donor [" + new DecimalFormat("0.##").format(fracRemoved * 100) + "%]");
  595. }
  597. return cellBarcodes;
  598. }
  600. /**
  601. * Parse the single donor likelihood file, retrieve the best donor for each cell barcode.
  602. *
  603. * @param singleDonorLikelihoodFile The OUTPUT file produced by AssignCellsToSamples.
  604. * @return A map where the key is the cell barcode, and the value is the best donor.
  605. */
  606. Map<String, String> getBestDonorForCell(final File singleDonorLikelihoodFile) {
  607. Map<String, String> result = new HashMap<>();
  609. CellCollectionSampleLikelihoodCollection cslc = CellCollectionSampleLikelihoodCollection.parseFile(singleDonorLikelihoodFile);
  611. for (String cellBarcode : cslc.getCellBarcodes()) {
  612. String donor = cslc.getLikelihoodCollection(cellBarcode).getBestSampleAssignment().getSample();
  613. result.put(cellBarcode, donor);
  614. }
  615. return (result);
  616. }
  618. // getNumberOfDonorsToTest(donorList.size(), this.MIN_DONOR_PAIRS_TESTED, this.FRACTION_DONOR_PAIRS_TESTED);
  619. private int getNumberOfDonorsToTest(int totalDonorPairs, int minDonorPairs, double fractionDonorPairs) {
  620. int numTest = (int) Math.ceil((double) totalDonorPairs * fractionDonorPairs);
  621. int result = Math.max(numTest, minDonorPairs);
  622. // don't test more than the total number of donors.
  623. // this is for when total donors are < minDonorPairs.
  624. if (result > totalDonorPairs)
  625. result = totalDonorPairs;
  626. log.info("Testing " + Integer.toString(result) + " Donor pairs per cell");
  627. return (result);
  629. }
  632. @Override
  633. protected String[] customCommandLineValidation() {
  634. IOUtil.assertFileIsReadable(this.VCF.toPath());
  635. this.INPUT_BAM = FileListParsingUtils.expandPicardHtsPathList(INPUT_BAM);
  636. IOUtil.assertFileIsWritable(this.OUTPUT);
  637. IOUtil.assertFileIsReadable(this.SAMPLE_FILE);
  638. IOUtil.assertFileIsReadable(this.SINGLE_DONOR_LIKELIHOOD_FILE);
  639. IOUtil.assertFileIsReadable(this.CELL_BC_FILE);
  640. final ArrayList<String> list = new ArrayList<>(1);
  641. if (OUTPUT_ALL_PAIRS != null)
  642. IOUtil.assertFileIsWritable(this.OUTPUT_ALL_PAIRS);
  643. if (OUTPUT_PER_SNP != null)
  644. IOUtil.assertFileIsWritable(this.OUTPUT_PER_SNP);
  646. if (!VCFUtils.hasIndex(this.VCF.toPath()))
  647. list.add("VCF is not indexed! Please index and retry.");
  649. if (this.CELL_CONTAMINATION_ESTIMATE_FILE != null && this.ALLELE_FREQUENCY_ESTIMATE_FILE == null)
  650. list.add("If CELL_CONTAMINATION_ESTIMATE_FILE is supplied, must also supply ALLELE_FREQUENCY_ESTIMATE_FILE");
  652. if (this.CELL_CONTAMINATION_ESTIMATE_FILE == null && this.ALLELE_FREQUENCY_ESTIMATE_FILE != null)
  653. list.add("If ALLELE_FREQUENCY_ESTIMATE_FILE is supplied, must also supply CELL_CONTAMINATION_ESTIMATE_FILE");
  655. if (this.FIXED_ERROR_RATE != null & this.MAX_ERROR_RATE == null)
  656. log.info("Running with a fixed error rate of " + this.FIXED_ERROR_RATE);
  658. if (this.MAX_ERROR_RATE != null && this.CELL_CONTAMINATION_ESTIMATE_FILE == null && this.ALLELE_FREQUENCY_ESTIMATE_FILE == null)
  659. log.info("Running with a maximum cap on error rate of " + this.MAX_ERROR_RATE);
  661. if (FRACTION_SAMPLES_PASSING < 0 | FRACTION_SAMPLES_PASSING > 1)
  662. list.add("FRACTION_SAMPLES_PASSING must be between 0 and 1, value was " + Double.toString(this.FRACTION_SAMPLES_PASSING));
  664. return CustomCommandLineValidationHelper.makeValue(super.customCommandLineValidation(), list);
  665. }
  667. /** Stock main method. */
  668. public static void main(final String[] args) {
  669. System.exit(new DetectDoublets().instanceMain(args));
  670. }
  672. /*
  673. private List<String> getExpectedSecondDonorsRankedByLikelihood(String cellBarcode, CellCollectionSampleLikelihoodCollection cslc, int pairsToTest,
  674. List<String> expectedDonors, String bestDonor) {
  675. // short circuit if the requested number of donors is the same as the number of expected donors IE the non-optimized
  676. // strategy.
  677. if (expectedDonors.size() == pairsToTest)
  678. return expectedDonors;
  680. CellSampleLikelihoodCollection c = cslc.getLikelihoodCollection(cellBarcode);
  681. if (c == null)
  682. throw new IllegalArgumentException("Could not find single donor likelihoods for cell " + cellBarcode);
  683. List<String> rankedDonors = c.getDonorsRankedByAssignmentLikelihood();
  684. // exclude the best donor from the ranked list. We don't want to select that.
  685. rankedDonors.remove(bestDonor);
  687. Set<String> expected = new HashSet<String>(expectedDonors);
  688. // filter ranked donors by expected.
  689. List<String> rankedExpectedDonors = rankedDonors.stream().filter(x -> expected.contains(x)).collect(Collectors.toList());
  691. // if pairs to test is set too high somehow (the sample lists doesn't match up with the input single donor assignments)
  692. // then limit the return.
  693. if (rankedExpectedDonors.size() < pairsToTest) {
  694. pairsToTest = rankedExpectedDonors.size();
  695. }
  696. // get the top <X> donors.
  697. rankedExpectedDonors = rankedExpectedDonors.subList(0, pairsToTest);
  698. return rankedExpectedDonors;
  700. }
  701. */
  704. }

DetectDoublets.java at commit d14776a, under MIT · at the source

Overview

Authors: Ricardo CH del Rosario1,2, Fenna M Krienen1,2, Qiangge Zhang2,3, Melissa Goldman1,2, Curtis Mello1,2, Alyssa Lutservitz1,2, Kiku Ichihara1,2, Alec Wysoker1,2, James Nemesh1,2, Guoping Feng2,3, Steven A McCarroll1,2,4
ORCID iDs: Ricardo CH del Rosario, Fenna M Krienen, Kiku Ichihara, Guoping Feng, Steven A McCarroll
  1. Department of Genetics, Harvard Medical School Boston United States
  2. Stanley Center for Psychiatric Research, Broad Institute of MIT and Harvard Cambridge United States
  3. McGovern Institute for Brain Research, Department of Brain and Cognitive Sciences, Massachusetts Institute of Technology Cambridge United States
  4. Howard Hughes Medical Institute Boston United States
Institutions: Harvard University (United States); Broad Institute (United States); Massachusetts Institute of Technology (United States); Howard Hughes Medical Institute (United States)
Journal: eLife, volume 13, article RP93640
Dates: published online 14 August 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.7554/elife.93640 · PMID 42598888 · PMCID PMC13476068 · OpenAlex W4393112346
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), non-human primate (organism), cellular / molecular (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Connectivity, Machine learning, fMRI & imaging
Keywords: Callithrix jacchus, sibling chimerism, microglia, marmosets, Other
MeSH: Callithrix*, Chimerism*, Microglia*, Animals, Brain, Female, Polymorphism, Single Nucleotide, Siblings (* major topic)
Journal subjects: Genetics and Genomics
Topic: Neuroinflammation and Neurodegeneration Mechanisms (Neurology, Neuroscience), according to OpenAlex
Funding: National Institutes of Health (U01MH114819); Stanley Center for Psychiatric Research, Broad Institute; James and Patricia Poitras Center for Psychiatric Disorders Research at MIT; Hock E. Tan and K. Lisa Yang Center for Autism Research at MIT
Citations: cited by 2 papers (Europe PMC); 34 references in the paper

Abstract

Chimerism happens rarely among most mammals, but is common in marmosets and tamarins, a result of fraternal twin or triplet birth patterns in which in utero connected circulatory systems (through which stem cells transit) lead to persistent blood chimerism (12–80%) throughout life. The presence of Y-chromosome DNA sequences in organs of female marmosets has long suggested that chimerism might also affect these organs. However, a longstanding question is whether this chimerism is driven by blood-derived cells or involves contributions from other cell types. To address this question, we analyzed single-cell RNA-seq data from blood, liver, kidney, and many brain regions across a number of marmosets, using transcribed single-nucleotide polymorphisms (SNPs) to identify cells with the sibling’s genome in various cell types within these tissues. Sibling-derived chimerism in all tissues arose entirely from cells of hematopoietic origin (i.e., myeloid and lymphoid lineages). In brain tissue this was reflected as sibling-derived chimerism among microglia (20–52%) and macrophages (18–64%) but not among other resident cell types (neurons, glia, or ependymal cells). The percentage of microglia that were sibling-derived showed significant variation across brain regions, even within individual animals, likely reflecting distinct responses by genetic-sibling microglia to local recruitment or proliferation cues or, potentially, distinct clonal expansion histories in different brain areas. In the animals and tissues we analyzed, microglial gene expression profiles bore a much stronger relationship to local/host context than to sibling genetic differences. Naturally occurring marmoset chimerism will provide new ways to recognize the effects of genes, mutations, and brain contexts on microglial biology and to distinguish between effects of microglia and other cell types on brain phenotypes.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 6 matches between paragraphs and lines of code.

broadinstitute/Drop-seq

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: d14776a599bbc401ef2d6ce5c6f96b23d0424cdf, 17 September 2026
Languages: Java (626), Python (58), R (19), Shell (14), Jupyter (1)
Size: 1,409 files, 718 scripts
Software Heritage: archived
Found in: “Software availability”
Holds: README, license file, environment (src/python/pyproject.toml, src/docker/java/Dockerfile, src/docker/PEER/Dockerfile, src/docker/python/Dockerfile, src/docker/R/Dockerfile, src/R/packages/DropSeq.dropulation/DESCRIPTION, src/R/packages/DropSeq.eqtl.susie/DESCRIPTION, src/R/packages/DropSeq.eqtl/DESCRIPTION, src/R/packages/DropSeq.utilities/DESCRIPTION), tests, continuous integration, documentation, 1 notebook
Not found: CITATION.cff
Tools: pandas (25 files), anndata (11 files), data.table (8 files), NumPy (8 files), SciPy (8 files), Scanpy (5 files), ggplot2 (3 files), Matplotlib (2 files), cowplot (1 file), scikit-learn (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
720 files
At the source: github.com/broadinstitute/Drop-seq

lh3/bwa

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: d82444c17edc2384420409f85557c6ae84019732, 7 August 2026
Languages: C (36), C/C++ (24), JavaScript (3), Perl (2), Shell (1)
Size: 80 files, 66 scripts
Software Heritage: archived
Found in: “Software availability”
Holds: README, license file, continuous integration
Not found: CITATION.cff, environment file, tests, documentation
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
68 files
At the source: github.com/lh3/bwa

satijalab/seurat

License: other
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 586015abde10618ecb32d3fe632267a83317a08d, 21 September 2026
Languages: R (114), C++ (8), C/C++ (4), C (1)
Size: 455 files, 127 scripts
Software Heritage: archived
Found in: “Software availability”
Holds: README, license file, environment (DESCRIPTION), tests, continuous integration, documentation, 70 notebooks
Not found: CITATION.cff
Tools: Seurat (75 files), ggplot2 (48 files), patchwork (27 files), tidyverse (19 files), cowplot (10 files), reshape2 (3 files), SingleCellExperiment (3 files), Plotly (2 files), data.table (1 file), DESeq2 (1 file), Harmony (1 file), igraph (1 file), limma (1 file), Monocle 3 (1 file), reticulate (1 file), UMAP (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
130 files
At the source: github.com/satijalab/seurat

IanevskiAleksandr/sc-type

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 630e15cf1e51f2612eda4ad0406dfb17503fa8c9, 10 November 2024
Languages: JavaScript (6), R (4)
Size: 117 files, 10 scripts
Software Heritage: not archived
Found in: “Software availability”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Seurat (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
12 files
At the source: github.com/IanevskiAleksandr/sc-type

satijalab/sctransform

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 49e35b5aeb76a602910207cbfde1561093340be3, 10 January 2026
Languages: R (28), C++ (2)
Size: 79 files, 30 scripts
Software Heritage: archived
Found in: “Software availability”
Holds: README, license file, environment (DESCRIPTION), tests, continuous integration, documentation, 11 notebooks
Not found: CITATION.cff
Tools: ggplot2 (12 files), reshape2 (11 files), tidyverse (8 files), Seurat (7 files), patchwork (4 files), cowplot (2 files), Plotly (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
32 files
At the source: github.com/satijalab/sctransform/

PMBio/peer

License: GPL-2.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 40bc4b2cd92459ce42f44dfe279717436395f3f6, 8 May 2012
Languages: C++ (869), C/C++ (648), Python (29), Shell (23), R (8), C (2)
Size: 2,526 files, 1,579 scripts
Software Heritage: not archived
Found in: “Software availability”
Holds: README, license file, environment (python/setup.py, R/peer/DESCRIPTION), documentation
Not found: CITATION.cff, tests, continuous integration
Tools: NumPy (2 files), SciPy (2 files), Matplotlib (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
1,582 files
At the source: github.com/PMBio/peer

Software availability

All software used in the analysis are publicly available. Drop-seq (analysis of snRNA-seq data, clustering, marker genes), Census-seq (estimation of chimerism in WGS data), and Dropulation analysis (estimation of chimerism in snRNA-seq data): https://github.com/broadinstitute/Drop-seq, Broad Institute, 2026; alignment and variant detection of Illumina WGS data: bwa (https://github.com/lh3/bwa, Li, 2026), GATK (https://gatk.broadinstitute.org), BCFtools (https://github.com/samtools/bcftools, Danecek, 2026), samtools (http://www.htslib.org/download), Picard Tools (https://broadinstitute.github.io/picard); R environment (https://www.rstudio.com/products/rstudio/download and https://www.r-project.org); single-cell analysis in R: Seurat (https://github.com/satijalab/seurat, Butler, 2026), cell type identification: scType (https://github.com/IanevskiAleksandr/sc-type, Aleksandr, 2024), SCT transform (https://github.com/satijalab/sctransform/, Hafemeister, 2026); PEER latent factor analysis (https://github.com/PMBio/peer, PMBio, 2012).

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 6 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 2,530 scripts, each with its path and the digest of its content;
  • 6 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

Brain snRNA-seq of 6 marmosets (CJ022, CJ023, CJ025, CJ026, CJ027, CJ028) were generated as part of the NIH's Brain Initiative Cell Census Network (BICCN) project, while brain snRNA-seq of 5 marmosets (CJ001, CJ006, CJ007, CJ023, CJ102), and all blood, liver, and kidney snRNA-seq were generated for this project. All snRNA-seq datasets are available in the BICCN NeMO portal (https://assets.nemoarchive.org/dat-hsgdsgu and https://assets.nemoarchive.org/dat-1je0mn3). The raw whole-genome sequencing datasets are available from the NIH Sequence Read Archive, under accession number BioProject PRJNA1068102.

The following datasets were generated:

del RosarioR 2023Marmosets Have their Birth Sibling's MicrogliaNeMOnemo:dat-hsgdsgu

Oregon Health and Science University 2024MCC Genome SequencingNCBI BioProjectPRJNA1068102

The following previously published datasets were used:

KrienenF 2023A marmoset brain cell census reveals persistent influence of developmental origin on neuronsNeMOnemo:dat-1je0mn3

del RosarioR 2023WGS of Marmoset Colony at Broad Institute by Stanley Center for Psychiatric ResearchNCBI Sequence Read ArchiveSRX32815946

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 11 authors, 5 keywords, 8 MeSH terms, 4 funders, 27 references.

Cite

This paper

del Rosario, R. C., Krienen, F. M., Zhang, Q., Goldman, M., Mello, C., Lutservitz, A., Ichihara, K., Wysoker, A., Nemesh, J., Feng, G., & McCarroll, S. A. (2026). Sibling chimerism among microglia in marmosets. eLife, 13, RP93640. https://doi.org/10.7554/elife.93640

BibTeX

@article{delrosario2026sibling,
author = {del Rosario, Ricardo CH and Krienen, Fenna M and Zhang, Qiangge and Goldman, Melissa and Mello, Curtis and Lutservitz, Alyssa and Ichihara, Kiku and Wysoker, Alec and Nemesh, James and Feng, Guoping and McCarroll, Steven A},
title = {{Sibling chimerism among microglia in marmosets}},
journal = {eLife},
year = {2026},
month = aug,
volume = {13},
pages = {RP93640},
publisher = {eLife Sciences Publications, Ltd},
issn = {2050-084X},
doi = {10.7554/elife.93640},
url = {https://doi.org/10.7554/elife.93640},
pmid = {42598888},
pmcid = {PMC13476068}
}

RIS

TY - JOUR
AU - del Rosario, Ricardo CH
AU - Krienen, Fenna M
AU - Zhang, Qiangge
AU - Goldman, Melissa
AU - Mello, Curtis
AU - Lutservitz, Alyssa
AU - Ichihara, Kiku
AU - Wysoker, Alec
AU - Nemesh, James
AU - Feng, Guoping
AU - McCarroll, Steven A
TI - Sibling chimerism among microglia in marmosets
T2 - eLife
J2 - Elife
PY - 2026
DA - 2026/08/14
VL - 13
SP - RP93640
SN - 2050-084X
PB - eLife Sciences Publications, Ltd
DO - 10.7554/elife.93640
UR - https://doi.org/10.7554/elife.93640
LA - en
ER -

CSL-JSON

{
"id": "10.7554/elife.93640",
"type": "article-journal",
"title": "Sibling chimerism among microglia in marmosets",
"container-title": "eLife",
"author": [
{
"family": "del Rosario",
"given": "Ricardo CH"
},
{
"family": "Krienen",
"given": "Fenna M"
},
{
"family": "Zhang",
"given": "Qiangge"
},
{
"family": "Goldman",
"given": "Melissa"
},
{
"family": "Mello",
"given": "Curtis"
},
{
"family": "Lutservitz",
"given": "Alyssa"
},
{
"family": "Ichihara",
"given": "Kiku"
},
{
"family": "Wysoker",
"given": "Alec"
},
{
"family": "Nemesh",
"given": "James"
},
{
"family": "Feng",
"given": "Guoping"
},
{
"family": "McCarroll",
"given": "Steven A"
}
],
"container-title-short": "Elife",
"volume": "13",
"page": "RP93640",
"DOI": "10.7554/elife.93640",
"PMID": "42598888",
"PMCID": "PMC13476068",
"ISSN": "2050-084X",
"publisher": "eLife Sciences Publications, Ltd",
"URL": "https://doi.org/10.7554/elife.93640",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
14
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s44318-026-00818-9 [code]
FAM134B-mediated ER-phagy degrades APP and suppresses Alzheimer's disease pathology.
Journal: The EMBO journal
In common: Monocle 3, Harmony, SingleCellExperiment, 20 other tools, cellular / molecular
[2] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Monocle 3, Harmony, SingleCellExperiment, 19 other tools, genetics / omics, cellular / molecular
[3] doi:10.1002/imt2.70163 [code]
Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: Monocle 3, Harmony, SingleCellExperiment, 18 other tools, genetics / omics, cellular / molecular
[4] doi:10.1016/j.xcrm.2026.102651 [code]
Integrative CSF profiling identifies disease-specific immune responses in leptomeningeal disease.
Journal: Cell reports. Medicine
In common: Harmony, SingleCellExperiment, reticulate, 18 other tools, genetics / omics, cellular / molecular
[5] doi:10.1038/s41586-026-10214-2 [code]
Multidimensional profiling of heterogeneity in supratentorial ependymomas.
Journal: Nature
In common: Monocle 3, Harmony, SingleCellExperiment, 17 other tools, genetics / omics
[6] doi:10.1002/ctm2.70683 [code]
Niacin promotes motor function recovery after spinal cord injury via Hcar2-dependent microglia immunometabolic regulation.
Journal: Clinical and translational medicine
In common: Monocle 3, Harmony, SingleCellExperiment, 13 other tools, non-human primate, genetics / omics, cellular / molecular
[7] doi:10.1186/s12974-026-03838-8 [code]
Acarbose modulates microglial Pkm2 acetylation to reshape immunometabolism and preserve retinal neurons after ischemia-reperfusion.
Journal: Journal of neuroinflammation
In common: Monocle 3, Harmony, SingleCellExperiment, 13 other tools, genetics / omics, cellular / molecular
[8] doi:10.1038/s41420-026-02971-w [code]
Multi-omics reveals heterogeneity and functional populations of oligodendrocyte progenitor cells induced by human neural stem cells.
Journal: Cell death discovery
In common: Monocle 3, Harmony, SingleCellExperiment, 13 other tools, genetics / omics, cellular / molecular
[9] doi:10.1038/s41467-026-73305-8 [code]
Comparative analysis of the cellular landscape in mammalian striatum.
Journal: Nature communications
In common: Harmony, SingleCellExperiment, anndata, 12 other tools, non-human primate, genetics / omics, cellular / molecular, 3 references
[10] doi:10.1038/s41467-026-76675-1 [code]
Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: Monocle 3, reticulate, limma, 16 other tools, genetics / omics

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.