OSCR

Single-cell profiling of DNA methylation in autism spectrum disorder prefrontal cortex reveals distinct regulatory and aging signatures.

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] § STAR★Methods › Quantification and statistical analysis › Bioinformatic alignment and feature quantification ↔ Notebooks/A06_mapping_star.ipynb, lines 16–63 · score 1.00 · BySJout, EndToEnd, alignEndsType, alignIntronMax, alignMatesGapMax, alignSJDBoverhangMin
  2. [2] § STAR★Methods › Quantification and statistical analysis › Bioinformatic alignment and feature quantification ↔ Notebooks/A06_mapping_star.ipynb, lines 545–630 · score 0.90 · pixel distance, optical duplicates, MarkDuplicates, samtools view, score, flags
  3. [3] § STAR★Methods › Quantification and statistical analysis › Bioinformatic alignment and feature quantification ↔ Notebooks/A04_mapping_bismark.ipynb, lines 425–501 · score 0.85 · pixel distance, optical duplicates, MarkDuplicates, conversion, samtools, mapped
  4. [4] § STAR★Methods › Quantification and statistical analysis › Bioinformatic alignment and feature quantification ↔ Notebooks/A03_trimming.ipynb, lines 57–196 · score 0.83 · adapter_fasta, cut_right, adapter sequences, snmCTseq, demultiplexing, trimmed
  5. [5] § STAR★Methods › Quantification and statistical analysis › Bioinformatic alignment and feature quantification ↔ Scripts/A02a_demultiplex_fastq.pl, lines 2–49 · score 0.53 · cell barcodes, demultiplexing, fasta, v0, nuclei, plates
  6. [6] § STAR★Methods › Quantification and statistical analysis › Enrichment analyses ↔ Notebooks/A00_environment_and_genome_setup.ipynb, lines 437–555 · score 0.51 · GENCODE v40, gtf, kb, downstream, exon, 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

Jupyter notebook · 1,120 lines · 33 KB · no license · 2 matches

  1. # %%
  2. # ## overall commands
  3. # ## might be cleaner to wait for STAR mapping to finish 100% (run A04a only --> check output)
  4. # ## in lieu of submitting all of the below, but technically could run all cmds at once
  6. # # * = job array based on "platenum"
  7. # # † = job array based on "batchnum" (two rows at a time)
  9. # qsub Scripts/A06a_star_mapping.sub # †
  10. # qsub Scripts/A06b_star_filtering.sub # †
  11. # qsub Scripts/A06c_check_star.sub
  12. # qsub Scripts/A06d_featurecounts.sub # *
  13. # qsub Scripts/A06e_star_bam_stats.sub # †
  16. # %% [markdown]
  17. # ## (A06a) star mapping
  19. # %%
  20. %%bash
  21. cat > ../Scripts/A06a_star_mapping.sub
  22. #!/bin/bash
  23. #$ -cwd
  24. #$ -o sublogs/A06a_star.$JOB_ID.$TASK_ID
  25. #$ -j y
  26. #$ -l h_rt=12:00:00,h_data=8G,exclusive
  27. #$ -pe shared 8
  28. #$ -N A06a_star
  29. #$ -t 1-512
  30. #$ -hold_jid_ad A03a_trim
  34. echo "Job $JOB_ID.$SGE_TASK_ID started on: " `hostname -s`
  35. echo "Job $JOB_ID.$SGE_TASK_ID started on: " `date `
  36. echo " "
  42. # environment init -------------------------------------------------------------
  44. . /u/local/Modules/default/init/modules.sh # <--
  45. module load anaconda3 # <--
  46. conda activate snmCTseq # <--
  48. export $(cat snmCT_parameters.env | grep -v '^#' | xargs) # <--
  50. skip_complete=true # <-- for help with incomplete jobs
  51. overwrite_partial=true # <-- for help with incomplete jobs
  53. # (most of the) STAR settings
  54. # originally adapted from ENCODE guidelines
  55. star_params="--runThreadN 8 \
  56. --genomeDir ${ref_starfolder} --genomeLoad LoadAndKeep \
  57. –alignEndsType EndToEnd --outSAMtype BAM Unsorted \
  58. --outSAMattributes NH HI AS NM MD --outSAMstrandField intronMotif \
  59. --sjdbOverhang 149 --outFilterType BySJout \
  60. --outFilterMultimapNmax 20 --alignSJoverhangMin 8 \
  61. --alignSJDBoverhangMin 1 –outFilterMismatchNmax 999 \
  62. --outFilterMismatchNoverLmax 0.04 --alignIntronMin 20 \
  63. --alignIntronMax 1000000 --alignMatesGapMax 1000000 --readFilesCommand zcat "
  67. # extract target filepaths -----------------------------------------------------
  69. # helper functions
  70. query_metadat () {
  71. awk -F',' -v targetcol="$1" \
  72. 'NR==1 {
  73. for (i=1;i<=NF;i++) {
  74. if ($i==targetcol) {assayout=i; break} }
  75. print $assayout
  76. }
  77. NR>1 {
  78. print $assayout
  79. }' ${metadat_well}
  80. }
  82. # extract target wells, print values for log
  83. batchnum=($(query_metadat "batchnum"))
  84. nwells=${#batchnum[@]}
  86. target_well_rows=()
  87. for ((row=1; row<=nwells; row++))
  88. do
  89. if [[ "${batchnum[$row]}" == "${SGE_TASK_ID}" ]]
  90. then
  91. target_well_rows+=($row)
  92. fi
  93. done
  97. # filepaths associated with target rows in well-level metadata -----------------
  99. wellprefix=($(query_metadat "wellprefix"))
  100. dir_well=($(query_metadat "A06a_dir_star"))
  102. # .fastqs for input to PE mapping (properly paired read pairs)
  103. fastq_r1p=($(query_metadat "A03a_fqgz_paired_R1"))
  104. fastq_r2p=($(query_metadat "A03a_fqgz_paired_R2"))
  106. # .fastqs for input to SE mapping, including singletons from trimming & unaligned in PE-mapping
  107. fastq_r1singletrim=($(query_metadat "A03a_fqgz_singletrim_R1"))
  108. fastq_r2singletrim=($(query_metadat "A03a_fqgz_singletrim_R2"))
  112. # temporary/intermediate mapping files -----------------------------------------
  113. # (expected output generated by STAR, given outname prefixes PE., SE1., SE2.)
  115. fastq_pe_unmap1=PE.Unmapped.out.mate1
  116. fastq_pe_unmap2=PE.Unmapped.out.mate2
  118. bam_pe=PE.Aligned.out.bam
  119. bam_se1=SE1.Aligned.out.bam
  120. bam_se2=SE2.Aligned.out.bam
  124. # run STAR mapping -------------------------------------------------------------
  126. if [[ ! -s mapping_star ]]
  127. then
  128. mkdir mapping_star
  129. fi
  131. cd ${dir_proj}/mapping_star
  133. # load genome index [5~10 min]
  134. # creates some apparent .log, .sam out despite just loading genome
  135. # so putting in mapping_star to keep these files in one place
  137. STAR --runThreadN 8 --genomeDir ${ref_starfolder} --genomeLoad LoadAndExit
  141. # loop through each well to map ------------------------------------------------
  143. for row in ${target_well_rows[@]}
  144. do
  146. # check directory/prior mapping --------------------------------------------
  147. cd ${dir_proj}
  149. # check for existing mapping output
  150. # if final outputs exist, skip; else run mapping .bam
  152. if [[ -s ${dir_well[$row]}/${bam_pe} \
  153. && -s ${dir_well[$row]}/${bam_se1} \
  154. && -s ${dir_well[$row]}/${bam_se2} \
  155. && "${skip_complete}"=="true" ]]
  156. then
  157. echo -e "final aligned .bams for '${wellprefix[$row]}' already exist. skipping this well.'"
  158. else
  159. echo -e "\n\napplying STAR to '${wellprefix[$row]}'...\n\n"
  161. # remove old directory if one exists to deal with incomplete files
  162. # albeit the only issue i've seen are .bai and .tbi indices
  163. # (these often are not overwritten by software in the pipeline,
  164. # resulting in "index is older than file" errors later on)
  165. if [[ -e ${dir_well[$row]} && "${overwrite_partial}" == "true" ]]
  166. then
  167. echo -e "\n\nWARNING: folder for '${wellprefix[$row]}' exists, but not its final .bam alignments."
  168. echo "because overwrite_partial=true, deleting the directory and re-mapping."
  169. rm -rf ${dir_well[$row]}
  170. fi
  172. mkdir ${dir_well[$row]}
  173. cd ${dir_well[$row]}
  177. # run alignments -----------------------------------------------------------
  178. # in: .fastqs from trimming: four .fastqs,
  179. # properly paired ($fastq_r2p, $fastq_r1p) and
  180. # trimming singletons ($fastq_r1singletrim, $fastq_r2singletrim)
  181. # out: - paired-end, single-end .bam alignments out (${bam_pe}, ${bam_se1}, ${bam_se2})
  182. # - key log files (e.g., mapping rate)
  183. # --------------------------------------------------------------------------
  185. # (i) paired-end mapping [<1-3 min] ----------------------------------------
  186. # assumptions: pairs that map ambiguously in paired-end mode should be discarded
  187. STAR ${star_params} \
  188. --outFileNamePrefix PE. \
  189. --readFilesIn ${dir_proj}/${fastq_r1p[$row]} ${dir_proj}/${fastq_r2p[$row]} \
  190. --outReadsUnmapped Fastx
  192. # .fq --> .fq.gz for future storage/help match STAR's expected input type [<1-2 min]
  193. bgzip ${fastq_pe_unmap1}
  194. bgzip ${fastq_pe_unmap2}
  198. # (ii.) single-end, R1 [<1-3 min] ------------------------------------------
  199. # includes Read 1 singletons from trimming and STAR mapping in (i)
  200. STAR ${star_params} \
  201. --outFileNamePrefix SE1. \
  202. --readFilesIn ${dir_proj}/${fastq_r1singletrim[$row]},${fastq_pe_unmap1}.gz
  206. # (iii.) single-end, R2 [<1-3 min] -----------------------------------------
  207. # includes Read 2 singletons from trimming and STAR mapping in (i)
  208. STAR ${star_params} \
  209. --outFileNamePrefix SE2. \
  210. --readFilesIn ${dir_proj}/${fastq_r2singletrim[$row]},${fastq_pe_unmap2}.gz
  214. # (iv.) optional clean-up (comment out as desired)
  215. # *.Log.final.out contains mapping rate; other *.out
  216. # have STAR internals that can be discarded
  218. rm -rf *_STARtmp
  219. rm -rf *SJ.out.tab
  220. rm "${fastq_pe_unmap1}.gz" "${fastq_pe_unmap2}.gz"
  223. fi
  224. done
  228. # unload genome from mem at end ------------------------------------------------
  230. STAR --genomeDir ${ref_starfolder} --genomeLoad Remove
  236. echo -e "\n\n'A06a_star_mapping' completed.\n\n"
  239. echo "Job $JOB_ID.$SGE_TASK_ID ended on: " `hostname -s`
  240. echo "Job $JOB_ID.$SGE_TASK_ID ended on: " `date `
  241. echo " "
  243. # %% [markdown]
  244. # ## (A06b) star filtering
  246. # %%
  247. %%bash
  248. cat > ../Scripts/A06b_classify_mCT_reads_STAR.pl
  250. #!/usr/bin/perl -w
  251. use strict;
  254. # A06b_classify_mCT_reads_STAR.pl, v0.2 =====================================================
  255. # based loosely on original perl script written by Dr. Chongyuan Luo (@luogenomics)
  256. # modifications by Choo Liu (@chooliu):
  257. # - readability/documentation
  258. # - changes in logic for paired-end mapping (fwd/rev strand) & more efficient parsing of MD:Z flag
  259. # incl looking at only C/G positions in read vs. ref changes, calculating # cytosines, etc
  260. # inputs: - .sam file from STAR (compatible with single-end or paired-end alignments)
  261. # outputs: - "_annotations" .tsv recording each alignment's:
  262. # number of cytosines, mC/C fraction, and call (DNA, RNA, ambiguous)
  263. # this annotation file is subsequently appended to the. bam to keep RNA reads
  264. # typical usage: perl A06b_classify_mCT_reads_STAR.pl alignments.sam
  265. # =====================================================================================
  267. # reading in .sam file
  268. my ($samfile)=($ARGV[0]);
  269. my @sample; my @samline;
  270. my $read; my @read; my @ref; my $ref; my @md; my $dir; my $unmch;
  272. # determining read classification
  273. my $filter_num_CH=3;
  274. my $filter_frac_mCHmin=0.5; # DNA (mCH/CH <0.5)
  275. my $filter_frac_mCHmax=0.9; # RNA (mCH/CH >0.9)
  276. my $mch_fraction; my $totalch; my $totcpg; my @call;
  281. # load .sam file, loop through lines in .sam
  282. # exports info on each read to "(input .sam file name)_annotations"
  284. # for each line in .sam file, --------------------------------------------------------
  285. # count unmethylated and methylated cytosines in CHG and CHH context based on XR-tag
  287. # notes:
  288. # - skip header header rows (starting with @)
  289. # - assumes the XR-tag is in column 10 (STAR output), $samline[9]
  290. # - currently ignores indels
  291. # - default settings: keep reads with >=3 CHNs and mCHN/CHN fraction >=0.9 [*]
  292. # at present, we only retain call="RNA" and discard all others
  294. open sam_in, "$samfile" or die $!;
  295. open sam_annotations, ">$samfile\_annotations" or die $!;
  296. while (<sam_in>)
  297. {
  298. chop $_;
  299. if (substr($_,0,1) eq '@') { print sam_annotations "\n"; }
  300. else
  301. {
  303. @samline = split(/\t/,$_);
  304. $read = $samline[9];
  305. @read = split(//,$read);
  306. @ref = @read;
  308. @md=split( /(\d+)/ , $samline[15]);
  310. # check mapping orientation relative to reference genome ----------------------------
  311. # dir=1 if fwd (99, 163, 0, 73), =0 if reverse (147, 83, 16, 89)
  312. # flaw in script is that these numbers manually specified--if see strange results,
  313. # check mapping output for other flags before filtering / add interpreter
  314. $dir = 0 + ( ($samline[1] eq 99) or ($samline[1] eq 163 ) or ($samline[1] eq 0) or ($samline[1] eq 73) );
  316. # compare read to reference genome --------------------------------------------------
  317. # split MD:Z: flag by numbers then examine resulting length
  318. # e.g., MD:Z:10A0B121 --> MD:Z: 10 A 0 B 121 has $num_md_features=5
  320. # assume fully methylated if perfect match to genomic ref
  321. # (if MD stores a single number indicating no changes, $#md = 1)
  322. my $unmch = 0;
  323. my $num_md_features=$#md;
  324. if ($num_md_features == 1) { }
  326. # <-- start ELSE for num_md_features != 1
  327. # otherwise, attempt to reconstruct reference sequence,
  328. # looping through positions in read where the read != ref nucleotide
  329. else {
  330. my $pos=0;
  332. for (my $i=1; $i<=$num_md_features; $i += 2) {
  334. # loop through bases where there are changes
  335. $pos = $pos + @md[$i];
  337. # do nothing if starts with "^" (deletion)
  338. # modify @ref sequence if differences
  339. if ( (rindex("@md[$i+1]", "^", 0) == 0) or ($i+1 > $num_md_features)) { }
  340. else {
  341. @ref[$pos] = @md[$i + 1];
  342. }
  344. # increment base by 1
  345. $pos += 1;
  346. }
  348. # tabulate # unmethylated cytosines ---------------------------------------------------
  349. # again looping through each position where read != ref
  350. # (although same loop as above, re-loop because potential cases where adjacent bases differ btwn ref & read)
  351. my $pos=0;
  352. for (my $i=1; $i<=$num_md_features; $i += 2) {
  354. $pos = $pos + @md[$i];
  356. # if read maps to forward strand of genome
  357. # check if at a CH-site, and unmethylated cytosine converted to "T"
  358. if ( ($dir eq 1) and
  359. (@ref[$pos] eq "C") and (@ref[$pos + 1] ne "G") and (@read[$pos] eq "T") ) {
  360. $unmch += 1;
  361. }
  363. # if read maps to reverse strand of genome
  364. # STAR rev compliments the $read sequence, so cytosines are represented by "G"
  365. # check if at a CH-site, and unmethylated cytosine converted to "A"
  366. # (note: mCH underestimated in niche case where dir=0 &
  367. # first base is cytosine, as @ref[$pos - 1] is undefined)
  368. if ( ($dir eq 0) and
  369. (@ref[$pos] eq "G") and (@ref[$pos - 1] ne "C") and (@read[$pos] eq "A") ) {
  370. $unmch += 1;
  371. }
  373. $pos += 1;
  374. }
  375. } # <-- end ELSE statement for num_md_features != 1
  377. $ref = join('', @ref);
  379. # count cytosines in CH-context --------------------------------------------------
  380. # note: since done via regex, faster to count [# C] and subtract [# CG] vs searching wild flag C[ACT]
  381. if ($dir eq 1) {
  382. my @totalch = $ref =~ /C/g;
  383. $totalch = scalar @totalch;
  384. my @totcpg = $ref =~ /CG/g;
  385. $totcpg = scalar @totcpg;
  386. }
  387. if ($dir eq 0) {
  388. my @totalch = $ref =~ /G/g;
  389. $totalch = scalar @totalch;
  390. my @totcpg = $ref =~ /CG/g;
  391. $totcpg = scalar @totcpg;
  392. }
  394. my $totalch = $totalch - $totcpg;
  396. # classify each read into modalities ----------------------------------------------
  397. if ($totalch==0) { # avoid div by zero error
  398. $mch_fraction=-999;
  399. @call="amb";
  400. } else {
  401. $mch_fraction = 1 - ($unmch/$totalch);
  402. if (($totalch>=$filter_num_CH) and ($mch_fraction < $filter_frac_mCHmin)) { @call="DNA"; }
  403. elsif (($totalch>=$filter_num_CH) and ($mch_fraction >= $filter_frac_mCHmax)) { @call="RNA"; }
  404. else { @call="amb"; } # exclude (low # CH or ambiguous mCH between 0.5-0.9)
  405. }
  407. print sam_annotations "${totalch}\t${mch_fraction}\t@call\n";
  409. }
  410. }
  412. close bam_in;
  414. # %%
  415. %%bash
  416. cat > ../Scripts/A06b_star_filtering.sub
  418. #!/bin/bash
  419. #$ -cwd
  420. #$ -o sublogs/A06b_starfilt.$JOB_ID.$TASK_ID
  421. #$ -j y
  422. #$ -l h_rt=6:00:00,h_data=16G
  423. #$ -N A06b_starfilt
  424. #$ -t 1-512
  425. #$ -hold_jid_ad A06a_star
  429. echo "Job $JOB_ID.$SGE_TASK_ID started on: " `hostname -s`
  430. echo "Job $JOB_ID.$SGE_TASK_ID started on: " `date `
  431. echo " "
  437. # environment init -------------------------------------------------------------
  439. . /u/local/Modules/default/init/modules.sh # <--
  440. module load anaconda3 # <--
  441. conda activate snmCTseq # <--
  443. export $(cat snmCT_parameters.env | grep -v '^#' | xargs) # <--
  445. skip_complete=true # <-- for help with incomplete jobs
  449. # extract target filepaths -----------------------------------------------------
  451. # helper functions
  452. query_metadat () {
  453. awk -F',' -v targetcol="$1" \
  454. 'NR==1 {
  455. for (i=1;i<=NF;i++) {
  456. if ($i==targetcol) {assayout=i; break} }
  457. print $assayout
  458. }
  459. NR>1 {
  460. print $assayout
  461. }' ${metadat_well}
  462. }
  464. # extract target wells, print values for log
  466. batchnum=($(query_metadat "batchnum"))
  467. nwells=${#batchnum[@]}
  469. target_well_rows=()
  470. for ((row=1; row<=nwells; row++))
  471. do
  472. if [[ "${batchnum[$row]}" == "${SGE_TASK_ID}" ]]
  473. then
  474. target_well_rows+=($row)
  475. fi
  476. done
  480. # filepaths associated with target rows in well-level metadata -----------------
  482. wellprefix=($(query_metadat "wellprefix"))
  483. dir_well=($(query_metadat "A06a_dir_star"))
  487. # set naming convention within each well folder --------------------------------
  489. # .bam files in (from A06a)
  490. bam_in_pe=PE.Aligned.out.bam
  491. bam_in_se1=SE1.Aligned.out.bam
  492. bam_in_se2=SE2.Aligned.out.bam
  494. # de-duplicated .bam
  495. bam_dedupe_pe=pe_dedupe.bam
  496. bam_dedupe_se1=se1_dedupe.bam
  497. bam_dedupe_se2=se2_dedupe.bam
  499. # de-duplication logs
  500. picard_log_pe=star_dedupe_pe.log
  501. picard_log_se1=star_dedupe_se1.log
  502. picard_log_se2=star_dedupe_se2.log
  504. # final .sam
  505. sam_q10_pe=q10_pe.sam
  506. sam_q10_se1=q10_se1.sam
  507. sam_q10_se2=q10_se2.sam
  509. # final .bam
  510. bam_final_pe=PE.Final.bam
  511. bam_final_se1=SE1.Final.bam
  512. bam_final_se2=SE2.Final.bam
  516. # map each well
  517. for row in ${target_well_rows[@]}
  518. do
  520. # check directory/prior mapping -------------------------------------------
  521. # check for existing mapping output
  522. # if final outputs exist, skip; else run filtering of direct STAR align
  524. cd ${dir_proj}
  526. if [[ -s ${dir_well[$row]}/${bam_final_pe} \
  527. && -s ${dir_well[$row]}/${bam_final_se1} \
  528. && -s ${dir_well[$row]}/${bam_final_se2} \
  529. && "${skip_complete}"=="true" ]]
  530. then
  531. echo -e "final aligned .bams for '${wellprefix[$row]}' already exist. skipping this well.'"
  532. else
  534. echo -e "\n\napplying STAR to '${wellprefix[$row]}'...\n\n"
  535. cd ${dir_proj}/${dir_well[$row]}
  538. # proceed if all input files exist
  539. if [[ -s ${bam_in_pe} && -s ${bam_in_se1} && -s ${bam_in_se2} ]]
  540. then
  542. # clear intermediate files (if not all of them exist)
  543. rm ${bam_final_pe} ${bam_final_se1} ${bam_final_se2} 2>&1 >/dev/null
  545. # run RNA filtering -------------------------------------------------------
  546. # in: three .bams from STAR mapping: $bam_in_* (paired-end, single-end read 1, read 2)
  547. # out: - MAPQ and f ($bam_final_*)
  548. # - key log files (e.g., mapping rate)
  549. # -------------------------------------------------------------------------
  551. # STAR will output singletons in .bam file
  552. # (e.g., read 2 maps but read 1 doesn't; samtools view -f 8 ${bam_in_pe})
  553. # hence the "-f 0x0002" flag for "proper pairs" only, and SE alignments merged with...
  554. # in later sections, "-f 0x0048" (the read is R1; it mapped but R2 mate didn't map)
  555. # and "-f 0x0088" (the read is R2; it mapped but R1 mate didn't map)
  557. # if both R1 & R2 of pair need to pass filtering criteria (AND instead of OR),
  558. # needs an added annotations --> wide step like the below
  559. # sed '$!N;s/\n/ /' ${sam_q10_pe}_annotations \
  560. # | awk '{print ( ($3 == "RNA") && ($6 == "RNA") )}' \
  561. # | awk '{print $0}1' > ${sam_q10_pe}_annotations_bothpairs
  562. # then awk filter off of this 'bothpairs' file
  564. # each step usually ~3-5 min
  566. # (i) paired-end
  567. picard MarkDuplicates -I ${bam_in_pe} \
  568. --ASSUME_SORT_ORDER "queryname" --OPTICAL_DUPLICATE_PIXEL_DISTANCE 2500 \
  569. --ADD_PG_TAG_TO_READS false --REMOVE_DUPLICATES \
  570. -O ${bam_dedupe_pe} -M ${picard_log_pe}
  572. samtools view -h -q 10 -f 0x0002 ${bam_dedupe_pe} > ${sam_q10_pe}
  573. perl ${dir_proj}/Scripts/A06b_classify_mCT_reads_STAR.pl ${sam_q10_pe}
  574. awk 'NR == FNR { if ($0=="" || $3=="RNA") line[NR]; next } (FNR in line)' \
  575. ${sam_q10_pe}_annotations ${sam_q10_pe} |
  576. samtools view -b - | samtools sort - > ${bam_final_pe}
  577. samtools index ${bam_final_pe}
  581. # (ii) single-end, read 1
  582. picard MarkDuplicates -I ${bam_in_se1} \
  583. --ASSUME_SORT_ORDER "queryname" --OPTICAL_DUPLICATE_PIXEL_DISTANCE 2500 \
  584. --ADD_PG_TAG_TO_READS false --REMOVE_DUPLICATES \
  585. -O ${bam_dedupe_se1} -M ${picard_log_se1}
  587. samtools view -h -q 10 ${bam_dedupe_se1} > ${sam_q10_se1}
  588. perl ${dir_proj}/Scripts/A06b_classify_mCT_reads_STAR.pl ${sam_q10_se1}
  589. awk 'NR == FNR { if ($0=="" || $3=="RNA") line[NR]; next } (FNR in line)' \
  590. ${sam_q10_se1}_annotations ${sam_q10_se1} |
  591. samtools view -b - | samtools sort - > ${bam_final_se1}
  592. samtools index ${bam_final_se1}
  596. # (iii) single-end, read 2
  597. picard MarkDuplicates -I ${bam_in_se2} \
  598. --ASSUME_SORT_ORDER "queryname" --OPTICAL_DUPLICATE_PIXEL_DISTANCE 2500 \
  599. --ADD_PG_TAG_TO_READS false --DUPLICATE_SCORING_STRATEGY RANDOM --REMOVE_DUPLICATES \
  600. -O ${bam_dedupe_se2} -M ${picard_log_se2}
  602. samtools view -h -q 10 ${bam_dedupe_se2} > ${sam_q10_se2}
  603. perl ${dir_proj}/Scripts/A06b_classify_mCT_reads_STAR.pl ${sam_q10_se2}
  604. awk 'NR == FNR { if ($0=="" || $3=="RNA") line[NR]; next } (FNR in line)' \
  605. ${sam_q10_se2}_annotations ${sam_q10_se2} |
  606. samtools view -b - | samtools sort - > ${bam_final_se2}
  607. samtools index ${bam_final_se2}
  611. # (iv.) optional clean-up (comment out as desired)
  612. rm ${bam_dedupe_pe} ${bam_dedupe_se1} ${bam_dedupe_se2}
  613. rm ${sam_q10_pe} ${sam_q10_se1} ${sam_q10_se2}
  614. # rm *annotations
  616. fi
  617. fi
  618. done
  624. echo -e "\n\n'A06b_star_filtering' completed.\n\n"
  628. echo "Job $JOB_ID.$SGE_TASK_ID ended on: " `hostname -s`
  629. echo "Job $JOB_ID.$SGE_TASK_ID ended on: " `date `
  630. echo " "
  632. # %% [markdown]
  633. # ## (A06b) check star output
  635. # %%
  636. %%bash
  637. cat > ../Scripts/A06c_check_star.sub
  639. #!/bin/bash
  640. #$ -cwd
  641. #$ -o sublogs/A06c_check_star.$JOB_ID
  642. #$ -j y
  643. #$ -l h_rt=1:00:00,h_data=8G
  644. #$ -N A06c_starcheck
  645. #$ -hold_jid A06a_star
  649. echo "Job $JOB_ID started on: " `hostname -s`
  650. echo "Job $JOB_ID started on: " `date `
  651. echo " "
  657. # environment init -------------------------------------------------------------
  659. export $(cat snmCT_parameters.env | grep -v '^#' | xargs) # <--
  663. # helper functions / extract target filepaths ----------------------------------
  665. query_metadat () {
  666. awk -F',' -v targetcol="$1" \
  667. 'NR==1 {
  668. for (i=1;i<=NF;i++) {
  669. if ($i==targetcol) {assayout=i; break} }
  670. print $assayout
  671. }
  672. NR>1 {
  673. print $assayout
  674. }' ${metadat_well}
  675. }
  677. check_filepaths_in_assay() {
  678. for file in $@
  679. do
  680. if [[ ! -s ${file} ]]
  681. then
  682. echo "missing '${file}'"
  683. fi
  684. done
  685. }
  687. check_filepath_by_batch() {
  688. target_array=($@)
  689. batches_to_rerun=()
  690. for ((target_batch=1; target_batch<=nbatches; target_batch++))
  691. do
  692. target_well_rows=()
  693. for ((row=1; row<=nwells; row++))
  694. do
  695. if [[ "${batchnum[$row]}" == "${target_batch}" ]]
  696. then
  697. target_well_rows+=($row)
  698. fi
  699. done
  701. batch_file_list=${target_array[@]: ${target_well_rows[0]}:${#target_well_rows[@]} }
  703. num_files_missing=$(check_filepaths_in_assay ${batch_file_list[@]} | wc -l)
  705. if [[ ${num_files_missing} > 0 ]]
  706. then
  707. batches_to_rerun+=(${target_batch})
  708. echo -e "${target_batch} \t ${num_files_missing}"
  709. fi
  710. done
  712. if [[ ${#batches_to_rerun[@]} > 0 ]]
  713. then
  714. echo "batches to re-run:"
  715. echo "${batches_to_rerun[*]}"
  716. fi
  717. }
  719. batchnum=($(query_metadat "batchnum"))
  721. nwells=${#batchnum[@]}
  722. nbatches=${batchnum[-1]}
  726. # apply checks for A06a output -------------------------------------------------
  728. echo "-----------------------------------------------------------------"
  729. echo "A. printing number of missing STAR mapping .bam files missing (by batch)... "
  730. echo "-----------------------------------------------------------------"
  732. bam_star_pe=($(query_metadat "A06a_bam_star_PE"))
  733. bam_star_se1=($(query_metadat "A06a_bam_star_SE1"))
  734. bam_star_se2=($(query_metadat "A06a_bam_star_SE2"))
  736. echo -e "\n\n\nchecking PE.Aligned.out.bam ----------------------------------------------------"
  737. echo -e "batchnum\tnum_missing"
  738. check_filepath_by_batch ${bam_star_pe[@]}
  740. echo -e "\n\n\nchecking SE1.Aligned.out.bam ---------------------------------------------------"
  741. echo -e "batchnum\tnum_missing"
  742. check_filepath_by_batch ${bam_star_se1[@]}
  744. echo -e "\n\n\nchecking SE2.Aligned.out.bam ---------------------------------------------------"
  745. echo -e "batchnum\tnum_missing"
  746. check_filepath_by_batch ${bam_star_se2[@]}
  748. echo "-----------------------------------------------------------------"
  749. echo "B. printing number of missing filtered .bam files missing (by batch)... "
  750. echo "-----------------------------------------------------------------"
  752. bam_starfilt_pe=($(query_metadat "A06b_bam_starfilt_PE"))
  753. bam_starfilt_se1=($(query_metadat "A06b_bam_starfilt_SE1"))
  754. bam_starfilt_se2=($(query_metadat "A06b_bam_starfilt_SE2"))
  756. echo -e "\n\n\nchecking PE.Final.out.bam ------------------------------------------------------"
  757. echo -e "batchnum\tnum_missing"
  758. check_filepath_by_batch ${bam_starfilt_pe[@]}
  760. echo -e "\n\n\nchecking SE1.Final.out.bam -----------------------------------------------------"
  761. echo -e "batchnum\tnum_missing"
  762. check_filepath_by_batch ${bam_starfilt_se1[@]}
  764. echo -e "\n\n\nchecking SE2.Final.out.bam -----------------------------------------------------"
  765. echo -e "batchnum\tnum_missing"
  766. check_filepath_by_batch ${bam_starfilt_se2[@]}
  768. echo -e "\n\n\nsuggest re-running and checking sublog output of above batches."
  772. echo -e "\n\n-----------------------------------------------------------------"
  773. echo "C. checking log files for issues."
  774. echo -e "-----------------------------------------------------------------\n"
  777. echo -e "\n\nchecking if 'completed' in sublogs/A06a_star* output."
  778. echo "if any filename is printed, the associated batch may have not completed mapping."
  780. grep -c "'A06a_star_mapping' completed" sublogs/A06a_star* | awk -F ":" '$2==0 {print $1}'
  783. echo -e "\n\nchecking if 'completed' is in sublogs/A06b_starfilt* output."
  784. echo "if any filename is printed, the associated batch may have not completed mapping."
  786. grep -c "A06b_star_filtering" sublogs/A06b_starfilt* | awk -F ":" '$2==0 {print $1}'
  789. echo -e "\n\nchecking if 'Exception' is in sublogs/A06b_starfilt* output."
  790. echo "if any filename is printed, may be issues with .bams from STAR mapping."
  791. echo "(e.g., stale file handle, premature end of file, unexpected compressed block length)"
  792. grep -c "Exception" sublogs/A06b_starfilt* | awk -F ":" '$2!=0 {print $1}'
  799. echo -e "\n\n'A06c_starcheck' completed.\n\n"
  802. echo "Job $JOB_ID ended on: " `hostname -s`
  803. echo "Job $JOB_ID ended on: " `date `
  804. echo " "
  806. # %% [markdown]
  807. # ## (A06d) feature counts
  809. # %%
  810. %%bash
  811. cat > ../Scripts/A06d_featurecounts.sub
  813. #!/bin/bash
  814. #$ -cwd
  815. #$ -o sublogs/A06d_featurecounts.$JOB_ID.$TASK_ID
  816. #$ -j y
  817. #$ -l h_rt=3:00:00,h_data=24G
  818. #$ -N A06d_featurecounts
  819. #$ -t 1-32
  820. #$ -hold_jid A06b_starfilt
  824. echo "Job $JOB_ID.$SGE_TASK_ID started on: " `hostname -s`
  825. echo "Job $JOB_ID.$SGE_TASK_ID started on: " `date `
  826. echo " "
  832. # environment init -------------------------------------------------------------
  834. . /u/local/Modules/default/init/modules.sh # <--
  835. module load anaconda3 # <--
  836. conda activate snmCTseq # <--
  838. export $(cat snmCT_parameters.env | grep -v '^#' | xargs) # <--
  840. dir_out_gene=featurecounts_gene/ # <--
  841. dir_out_exon=featurecounts_exon/ # <--
  843. quantify_exons="true" # <-- (typically we analyze combined exon+intron)
  847. # extract target filepaths -----------------------------------------------------
  849. # helper functions
  850. query_metadat () {
  851. awk -F',' -v targetcol="$1" \
  852. 'NR==1 {
  853. for (i=1;i<=NF;i++) {
  854. if ($i==targetcol) {assayout=i; break} }
  855. print $assayout
  856. }
  857. NR>1 {f
  858. print $assayout
  859. }' ${metadat_well}
  860. }
  862. # extract target wells, print values for log
  864. platenum=($(query_metadat "platenum"))
  865. nwells=${#platenum[@]}
  867. target_well_rows=()
  868. for ((row=1; row<=nwells; row++))
  869. do
  870. if [[ "${platenum[$row]}" == "${SGE_TASK_ID}" ]]
  871. then
  872. target_well_rows+=($row)
  873. fi
  874. done
  878. # filepaths associated with target rows in well-level metadata -----------------
  880. wellprefix=($(query_metadat "wellprefix"))
  881. dir_well=($(query_metadat "A06a_dir_star"))
  883. bam_pe=($(query_metadat "A06b_bam_starfilt_PE"))
  884. bam_se1=($(query_metadat "A06b_bam_starfilt_SE1"))
  885. bam_se2=($(query_metadat "A06b_bam_starfilt_SE2"))
  889. # extract valid .bam filepaths -------------------------------------------------
  891. # checks only if the .bam exists (otherwise featureCounts may terminate with error)
  892. # could also wrap in check that number of alignments in file is > 0
  893. # for greater future compatibility with featureCounts
  895. pe_files=$(
  896. for file in ${bam_pe[@]: ${target_well_rows[0] }:${#target_well_rows[@] } }
  897. do
  898. if [[ -s ${file} ]]
  899. then
  900. echo ${file}
  901. fi
  902. done)
  904. se1_files=$(
  905. for file in ${bam_se1[@]: ${target_well_rows[0] }:${#target_well_rows[@] } }
  906. do
  907. if [[ -s ${file} ]]
  908. then
  909. echo ${file}
  910. fi
  911. done)
  913. se2_files=$(
  914. for file in ${bam_se2[@]: ${target_well_rows[0] }:${#target_well_rows[@] } }
  915. do
  916. if [[ -s ${file} ]]
  917. then
  918. echo ${file}
  919. fi
  920. done)
  924. # featurecounts on genes -------------------------------------------------------
  925. # usually <1 min/well --> <1 hr/plate
  927. if [[ ! -s ${dir_out_gene} ]]
  928. then
  929. mkdir ${dir_out_gene}
  930. fi
  932. echo "running gene featureCounts on paired-end alignments (mapping_star/*/PE.Final.bam)."
  933. featureCounts -p -T 4 -t gene -a ${ref_gtf} \
  934. -o ${dir_out_gene}/PE_${SGE_TASK_ID} --donotsort --tmpDir ${dir_scratch} ${pe_files}
  936. echo "running gene featureCounts on single-end R1 alignments (mapping_star/*/SE1.Final.bam)."
  937. featureCounts -T 4 -t gene -a ${ref_gtf} \
  938. -o ${dir_out_gene}/SE1_${SGE_TASK_ID} --donotsort --tmpDir ${dir_scratch} ${se1_files}
  940. echo "running gene featureCounts on single-end R2 alignments (mapping_star/*/SE2.Final.bam)."
  941. featureCounts -T 4 -t gene -a ${ref_gtf} \
  942. -o ${dir_out_gene}/SE2_${SGE_TASK_ID} --donotsort --tmpDir ${dir_scratch} ${se2_files}
  946. # featurecounts on exons -------------------------------------------------------
  947. # usually <1 min/well --> <1 hr/plate
  949. if [[ ${quantify_exons} == "true" ]]
  950. then
  952. if [[ ! -s ${dir_out_exon} ]]
  953. then
  954. mkdir ${dir_out_exon}
  955. fi
  957. echo "running exon featureCounts on paired-end alignments (mapping_star/*/PE.Final.bam)."
  958. featureCounts -p -T 4 -t exon -a ${ref_gtf} \
  959. -o ${dir_out_exon}/PE_${SGE_TASK_ID} --donotsort --tmpDir ${dir_scratch} ${pe_files}
  961. echo "running exon featureCounts on single-end R1 alignments (mapping_star/*/SE1.Final.bam)."
  962. featureCounts -T 4 -t exon -a ${ref_gtf} \
  963. -o ${dir_out_exon}/SE1_${SGE_TASK_ID} --donotsort --tmpDir ${dir_scratch} ${se1_files}
  965. echo "running exon featureCounts on single-end R2 alignments (mapping_star/*/SE2.Final.bam)."
  966. featureCounts -T 4 -t exon -a ${ref_gtf} \
  967. -o ${dir_out_exon}/SE2_${SGE_TASK_ID} --donotsort --tmpDir ${dir_scratch} ${se2_files}
  969. fi
  975. echo -e "\n\n'A06d_featurecounts' completed.\n\n"
  979. echo "Job $JOB_ID.$SGE_TASK_ID ended on: " `hostname -s`
  980. echo "Job $JOB_ID.$SGE_TASK_ID ended on: " `date `
  981. echo " "
  985. # %% [markdown]
  986. # ## (A06e) final .bam stats like counts, % intergenic
  988. # %%
  989. %%bash
  990. cat > ../Scripts/A06e_star_bam_stats.sub
  992. #!/bin/bash
  993. #$ -cwd
  994. #$ -o sublogs/A06e_samstat_star.$JOB_ID.$TASK_ID
  995. #$ -j y
  996. #$ -l h_rt=2:00:00,h_data=24G
  997. #$ -N A06e_samstat_star
  998. #$ -t 1-512
  999. #$ -hold_jid_ad A06b_starfilt
  1003. echo "Job $JOB_ID.$SGE_TASK_ID started on: " `hostname -s`
  1004. echo "Job $JOB_ID.$SGE_TASK_ID started on: " `date `
  1005. echo " "
  1011. # environment init -------------------------------------------------------------
  1013. . /u/local/Modules/default/init/modules.sh # <--
  1014. module load anaconda3 # <--
  1015. conda activate snmCTseq # <--
  1017. export $(cat snmCT_parameters.env | grep -v '^#' | xargs) # <--
  1019. skip_complete=true # <-- for help with incomplete jobs
  1023. # extract target filepaths -----------------------------------------------------
  1025. # helper functions
  1026. query_metadat () {
  1027. awk -F',' -v targetcol="$1" \
  1028. 'NR==1 {
  1029. for (i=1;i<=NF;i++) {
  1030. if ($i==targetcol) {assayout=i; break} }
  1031. print $assayout
  1032. }
  1033. NR>1 {
  1034. print $assayout
  1035. }' ${metadat_well}
  1036. }
  1038. # extract target wells, print values for log
  1040. batchnum=($(query_metadat "batchnum"))
  1041. nwells=${#batchnum[@]}
  1043. target_well_rows=()
  1044. for ((row=1; row<=$nwells; row++))
  1045. do
  1046. if [[ "${batchnum[$row]}" == "${SGE_TASK_ID}" ]]
  1047. then
  1048. target_well_rows+=($row)
  1049. fi
  1050. done
  1054. # filepaths associated with target rows in well-level metadata -----------------
  1056. wellprefix=($(query_metadat "wellprefix"))
  1057. dir_well=($(query_metadat "A06a_dir_star"))
  1061. # samtools stats on each well in the batch -------------------------------------
  1062. # usually <5 seconds per .bam --> <12 min/24 wells
  1064. for row in ${target_well_rows[@]}
  1065. do
  1067. cd ${dir_proj}
  1069. if [[ -s ${dir_well[$row]}/samstats_PE \
  1070. && -s ${dir_well[$row]}/samstats_SE1 \
  1071. && -s ${dir_well[$row]}/samstats_SE2 \
  1072. && -s ${dir_well[$row]}/picard_PE \
  1073. && -s ${dir_well[$row]}/picard_SE1 \
  1074. && -s ${dir_well[$row]}/picard_SE2 ]]
  1075. then
  1076. echo -e "final metrics for '${wellprefix[$row]}' already exist."
  1077. if [[ ${skip_complete} = "true" ]]
  1078. then
  1079. echo "skip_complete = true. skipping this well."
  1080. continue
  1081. else
  1082. echo "but skip_complete = false. re-running this well."
  1083. fi
  1084. fi
  1086. echo -e "\n\ngetting metrics for '${wellprefix[$row]}'...\n\n"
  1088. cd ${dir_well[$row]}
  1090. # run samtools stats
  1091. samtools stats PE.Final.bam | grep '^SN' | cut -f 2,3 > samstats_PE
  1092. samtools stats SE1.Final.bam | grep '^SN' | cut -f 2,3 > samstats_SE1
  1093. samtools stats SE2.Final.bam | grep '^SN' | cut -f 2,3 > samstats_SE2
  1095. # run picard CollectRnaSeqMetrics
  1096. samtools sort PE.Final.bam | \
  1097. picard CollectRnaSeqMetrics -I /dev/stdin -O picard_PE \
  1098. --REF_FLAT ${ref_flat} -STRAND "NONE" --RIBOSOMAL_INTERVALS ${ref_rrna}
  1099. samtools sort SE1.Final.bam | \
  1100. picard CollectRnaSeqMetrics -I /dev/stdin -O picard_SE1 \
  1101. --REF_FLAT ${ref_flat} -STRAND "NONE" --RIBOSOMAL_INTERVALS ${ref_rrna}
  1102. samtools sort SE2.Final.bam | \
  1103. picard CollectRnaSeqMetrics -I /dev/stdin -O picard_SE2 \
  1104. --REF_FLAT ${ref_flat} -STRAND "NONE" --RIBOSOMAL_INTERVALS ${ref_rrna}
  1106. done
  1112. echo -e "\n\n'A06e_samstat_star' completed.\n\n"
  1116. echo "Job $JOB_ID.$SGE_TASK_ID ended on: " `hostname -s`
  1117. echo "Job $JOB_ID.$SGE_TASK_ID ended on: " `date `
  1118. echo " "

A06_mapping_star.ipynb at commit 3a5f663, no license · at the source

Overview

Authors: Katherine W Eyring1, Cuining Liu2,3, Nasser Elhajjaoui3, Kevin D Abuhanna3, Yi Zhang3, Zachary von Behren3, Min Jen Tsai3, Eleazar Eskin3,4,5, Daniel H Geschwind1,3,6,7,8, Chongyuan Luo2,3
ORCID iDs: Zachary von Behren
  1. Department of Neurology, David Geffen School of Medicine, University of California, Los Angeles, Los Angeles, CA 90095, USA
  2. Bioinformatics Interdepartmental Program, University of California, Los Angeles, Los Angeles, CA, USA
  3. Department of Human Genetics, David Geffen School of Medicine, University of California, Los Angeles, Los Angeles, CA 90095, USA
  4. Department of Computer Science, University of California, Los Angeles, Los Angeles, CA 90095, USA
  5. Department of Computational Medicine, University of California, Los Angeles, Los Angeles, CA 90095, USA
  6. Center for Autism Research and Treatment, Semel Institute, David Geffen School of Medicine, University of California, Los Angeles, Los Angeles, CA 90095, USA
  7. Department of Psychiatry, Semel Institute, David Geffen School of Medicine, University of California, Los Angeles, Los Angeles, CA 90095, USA
  8. Institute for Precision Health, University of California, Los Angeles, Los Angeles, CA 90095, USA
Institutions: University of California, Los Angeles (United States)
Journal: Cell genomics, volume 6, issue 7, article 101278
Dates: received 27 July 2025; accepted 19 May 2026; published online 19 June 2026; in print July 2026
Type: Research article · Language: English
License: CC BY-NC
Identifiers: DOI 10.1016/j.xgen.2026.101278 · PMID 42320469 · PMCID PMC13347942 · OpenAlex W7165370473
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism), autism (population)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning
Keywords: neuroscience, neurodevelopment, autism, ASD, DNA methylation, single-cell genomics, multi-omics
MeSH: Aging*, Autism Spectrum Disorder*, DNA Methylation*, Prefrontal Cortex*, Single-Cell Analysis*, Adolescent, Adult, Child, Child, Preschool, Epigenesis, Genetic, Female, Humans, Male, Neurons, Young Adult (* major topic)
Topic: Autism Spectrum Disorder Research (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: NIMH NIH HHS (R01 MH125252, U01 MH130995); National Institute of Mental Health
Citations: cited by 1 paper (Europe PMC); 139 references in the paper
Research resources: Illumina Novaseq 6000 RRID:SCR_016387, Beckman Biomek i7 liquid handler RRID:SCR_018094, TC20 Automated Cell Counter RRID:SCR_025462

Abstract

Autism spectrum disorder (ASD) is a heterogeneous neurodevelopmental condition. Studies of postmortem ASD brain tissue have revealed convergent molecular changes across the cortex. Whether these features are reflected in cell-type-specific epigenetic signatures is unknown. Here, we present a single-cell analysis of DNA methylation (DNAm) coupled with transcriptomics in ASD. Using snmCT-seq, we profiled DNAm and transcript levels from over 60,000 nuclei derived from the prefrontal cortex of 49 donors. We identified over 30,000 differentially methylated regions (DMRs) in ASD that were enriched in promoters and cell-type-specific regulatory elements active across the lifespan. ASD-related methylation changes were uncorrelated with transcript levels and were small in magnitude compared with age-associated effects. Age-DMRs were concentrated in excitatory neurons and revealed distinct roles for CG and non-CG methylation. Age-varying methylation signatures of ASD identified neuron projection development as a key process perturbed in ASD, highlighting the heterogeneous impact of ASD across the lifespan.

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

Repository

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

chooliu/snmCTseq_Pipeline

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 3a5f663678849c9cdd1f666d79135b11e9c19b70, 17 March 2024
Languages: Python (18), Jupyter (9), Perl (3)
Size: 69 files, 30 scripts
Software Heritage: not archived
Found in: “Data and code availability”
Holds: README, documentation, 9 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration
Tools: pandas (18 files), NumPy (4 files), SAMtools (3 files), FastQC (2 files), STAR (2 files), Subread (featureCounts) (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
31 files
At the source: github.com/chooliu/snmCTseq_Pipeline

The paper's code and data availability statement is in the Data section.

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:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 30 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

Datasets cited

Data and code availability

All original code has been deposited at https://github.com/chooliu/snmCTseq_Pipeline and is publicly available as of the date of publication.

Processed data have been deposited at GEO as series no. GSE298243 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE298243) and are publicly available as of June 17, 2025.

Controlled-access data are available here: https://assets.nemoarchive.org/dat-htywcnb.

Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Reproduced under the paper's license (CC BY-NC), 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, issue, pages, dates, 10 authors, 7 keywords, 15 MeSH terms, 2 funders, 130 references, 3 RRIDs.

Cite

This paper

Eyring, K. W., Liu, C., Elhajjaoui, N., Abuhanna, K. D., Zhang, Y., von Behren, Z., Tsai, M. J., Eskin, E., Geschwind, D. H., & Luo, C. (2026). Single-cell profiling of DNA methylation in autism spectrum disorder prefrontal cortex reveals distinct regulatory and aging signatures. Cell genomics, 6(7), 101278. https://doi.org/10.1016/j.xgen.2026.101278

BibTeX

@article{eyring2026single,
author = {Eyring, Katherine W and Liu, Cuining and Elhajjaoui, Nasser and Abuhanna, Kevin D and Zhang, Yi and von Behren, Zachary and Tsai, Min Jen and Eskin, Eleazar and Geschwind, Daniel H and Luo, Chongyuan},
title = {{Single-cell profiling of DNA methylation in autism spectrum disorder prefrontal cortex reveals distinct regulatory and aging signatures}},
journal = {Cell genomics},
year = {2026},
month = jun,
volume = {6},
number = {7},
pages = {101278},
publisher = {Elsevier},
issn = {2666-979X},
doi = {10.1016/j.xgen.2026.101278},
url = {https://doi.org/10.1016/j.xgen.2026.101278},
pmid = {42320469},
pmcid = {PMC13347942}
}

RIS

TY - JOUR
AU - Eyring, Katherine W
AU - Liu, Cuining
AU - Elhajjaoui, Nasser
AU - Abuhanna, Kevin D
AU - Zhang, Yi
AU - von Behren, Zachary
AU - Tsai, Min Jen
AU - Eskin, Eleazar
AU - Geschwind, Daniel H
AU - Luo, Chongyuan
TI - Single-cell profiling of DNA methylation in autism spectrum disorder prefrontal cortex reveals distinct regulatory and aging signatures
T2 - Cell genomics
J2 - Cell Genom
PY - 2026
DA - 2026/06/19
VL - 6
IS - 7
SP - 101278
SN - 2666-979X
PB - Elsevier
DO - 10.1016/j.xgen.2026.101278
UR - https://doi.org/10.1016/j.xgen.2026.101278
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.xgen.2026.101278",
"type": "article-journal",
"title": "Single-cell profiling of DNA methylation in autism spectrum disorder prefrontal cortex reveals distinct regulatory and aging signatures",
"container-title": "Cell genomics",
"author": [
{
"family": "Eyring",
"given": "Katherine W"
},
{
"family": "Liu",
"given": "Cuining"
},
{
"family": "Elhajjaoui",
"given": "Nasser"
},
{
"family": "Abuhanna",
"given": "Kevin D"
},
{
"family": "Zhang",
"given": "Yi"
},
{
"family": "von Behren",
"given": "Zachary"
},
{
"family": "Tsai",
"given": "Min Jen"
},
{
"family": "Eskin",
"given": "Eleazar"
},
{
"family": "Geschwind",
"given": "Daniel H"
},
{
"family": "Luo",
"given": "Chongyuan"
}
],
"container-title-short": "Cell Genom",
"volume": "6",
"issue": "7",
"page": "101278",
"DOI": "10.1016/j.xgen.2026.101278",
"PMID": "42320469",
"PMCID": "PMC13347942",
"ISSN": "2666-979X",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.xgen.2026.101278",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
19
]
]
}
}

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/s41593-026-02247-7 [code]
Transcriptomic and phenotypic convergence of neurodevelopmental disorder risk genes in vitro and in vivo.
Journal: Nature neuroscience
In common: autism, genetics / omics, 12 references
[2] doi:10.1038/s42003-026-10059-5 [code]
Transcriptomic analysis in autism spectrum disorder suggests three molecular subtypes with distinct phenotypic profiles and functional pathways.
Journal: Communications biology
In common: FastQC, STAR, SAMtools, autism, genetics / omics, 6 references
[3] doi:10.1016/j.celrep.2026.117073 [code]
Single-cell epigenomics uncovers heterochromatin instability and transcription factor dysfunction during mouse brain aging.
Journal: Cell reports
In common: Subread (featureCounts), SAMtools, pandas, 1 other tool, genetics / omics, 7 references
[4] doi:10.1016/j.xhgg.2026.100652 [code]
CRISPR-engineered deletion of POGZ alters transcription factor binding at promoters of genes involved in synaptic signaling.
Journal: HGG advances
In common: STAR, SAMtools, pandas, 1 other tool, autism, 8 references
[5] 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: STAR, SAMtools, pandas, 1 other tool, autism, genetics / omics, 7 references
[6] doi:10.21203/rs.3.rs-9927928/v1 [code]
Genome-wide and allele-resolved maps of the radial architecture of the mouse genome
Journal: Research Square (preprint)
In common: FastQC, Subread (featureCounts), STAR, 3 other tools, genetics / omics, 4 references
[7] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: FastQC, Subread (featureCounts), STAR, 3 other tools, genetics / omics, 4 references
[8] doi:10.1038/s41467-026-74320-5 [code]
Spatial architecture of autism pathogenesis reveals mosaic structural disarray during early development.
Journal: Nature communications
In common: pandas, NumPy, autism, genetics / omics, 9 references
[9] doi:10.1038/s41586-026-10512-9 [code]
Astrocyte glucocorticoid receptor signalling restricts neuronal plasticity.
Journal: Nature
In common: Subread (featureCounts), STAR, SAMtools, 2 other tools, 6 references
[10] doi:10.1038/s41586-026-10679-1 [code]
Cortical development dynamics across autism spectrum disorder mouse models.
Journal: Nature
In common: pandas, NumPy, autism, genetics / omics, 10 references

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.