Update ATAC-seq pipeline by marinafloresp · Pull Request #322 · ejh243/BrainFANS

marinafloresp · 2025-09-05T14:32:33Z

Description

This pull request will update the ATAC-seq pipeline with the latest working scripts. Documentation still needs to be fully updated. The scripts of this update still need to be reviewed.

Issue ticket number

This pull request is to address issue: #number.

Type of pull request

Bug fix
New feature/enhancement
Code refactor
Documentation update

Checklist

I have performed a self-review of my own code
I have commented my code, particularly in hard-to-understand areas
I have tested my code to check that it is functional
I have used linters to check for common sources of errors
I have implemented fail safes in my code to account for edge cases
I have made the corresponding changes to the documentation

sof202

Some nice changes here Marina. Some of the explanations in the Rmarkdown have much more clarity and succinctness than before.

This was pretty hard to review however. All I see on my end is one lump sum of all of your changes from the last couple of months. It would be easier to see your intentions behind changes if each change came as a different commit. As such I must assume I have missed a lot of minor details here. For example, you might remove a function somehwere, but I have no idea why it is has been removed (perhaps it has been replaced, moved or simply isn't needed anymore, I have no way of knowing).

I trust that you and Philippa have already been using these scripts and so they should be functional.

sof202 · 2025-09-05T14:40:55Z

 metadata<-read.table(sampleSheet, header = TRUE, sep = ',', stringsAsFactors = FALSE)
-
-samplesPassed <- read.csv(file.path(paste0(metaDir,"/passAllStatus.csv")), sep = ",",stringsAsFactors = FALSE)
-samplesPassed <- samplesPassed[samplesPassed$QCS1 == TRUE & samplesPassed$QCS2 == TRUE,]
-qcStats <- read.csv(file.path(paste0(metaDir,"/samplesGroupAnalysisQC1.csv")), sep = ",",stringsAsFactors = FALSE)
-
-if(file.exists(paste0(metaDir,"/leaveOutSamples.txt"))){
-  loSamples <- read.table(file.path(paste0(metaDir,"/leaveOutSamples.txt")))[,1]
-  samplesStatus <- samplesPassed[!samplesPassed$sampleID %in% loSamples,]
-  pheno <- metadata[metadata$sampleID %in% samplesStatus$sampleID,]
-  samplesStatus <- samplesStatus[samplesStatus$sampleID %in% pheno$sampleID,]
-  qcStats <- qcStats[qcStats$Sample %in% samplesStatus$sampleID,]
-}else{
-  pheno <- metadata[metadata$sampleID %in% samplesPassed$sampleID,]
-  qcStats <- qcStats[qcStats$Sample %in% samplesPassed$sampleID,]
-
-}
-
-df_ordered <- pheno[order(pheno$fraction), ]
+samplesPassed <- read.csv(file.path(qc12stats), sep = ",",stringsAsFactors = FALSE)
+samples <- read.csv(file.path(groupAnalysisSamples), sep = ",",stringsAsFactors = FALSE)
+samplesPassed <- samplesPassed[samplesPassed$PASSALL == TRUE,]
+qcStats <- read.table(qc1stats, header = TRUE, sep = ',', stringsAsFactors = FALSE)


Can I ask why all of our scrips state stringsAsFactors = FALSE. Do we ever have this as TRUE anywhere at all? The default behaviour is this being set to FALSE. The separator being set here is also not really required as it is also the default, but I'm more ok with this as it make it more obvious what the expected format is.

sof202 · 2025-09-05T14:42:40Z

-}
-
-df_ordered <- pheno[order(pheno$fraction), ]
+samplesPassed <- read.csv(file.path(qc12stats), sep = ",",stringsAsFactors = FALSE)


Is file.path() required here? I assume that qc12stats is just some string, in which case this function call is not doing anything. If it qc12stats is a vector of strings, then nevermind.

sof202 · 2025-09-05T14:51:59Z

+samplesPassed <- samplesPassed[samplesPassed$PASSALL == TRUE,]
+qcStats <- read.table(qc1stats, header = TRUE, sep = ',', stringsAsFactors = FALSE)
+metadata$sequencingBatch <- as.factor(metadata$sequencingBatch)
+metadata<- droplevels(metadata)


What am I missing here? Is this actually doing anything? This call will remove any levels in a factor vector that has no representation in the vector:

a <- factor(c(1,1,2,2,3)) a[5] <- 2 # 3 is still a level despite not appearing in a a <- droplevels(a) # 3 is removed as a level

But here, the only factor column is the sequencingBatch column. The levels have just been generated in the previous line and so all levels appear in the column still. This therefore does nothing right?

sof202 · 2025-09-05T14:53:14Z

@@ -81,8 +67,6 @@ colors <- brewer.pal(3, "Set2")
 ## Load counts in peaks files and match to pheno samples data


stale comment?

sof202 · 2025-09-05T14:58:31Z

  tpm.scores<-tpm.scores[,which(tpm.pca > 0.01)]
  percentVar = round(100 * (tpm.pca) )
-  results <- list(tpm.scores[,c(1:5)], percentVar)
+  results <- list(tpm.scores[,c(1:4)], percentVar)


Any particular reason we are now taking 4 instead of 5 principle components? If this is something that would end up being variable, it might be a good idea to add it as a function parameter so the caller can decide how many to pick.

sof202 · 2025-09-05T16:02:50Z

+## - <cell-group>.filt.broadPeak.bed" in 4_calledpeaks folder                                                     ||
+##                                                                                                                ||
+## REQUIRES:                                                                                                      ||
+## - R version > 4.3                                                                                              ||


Here we now require a version of R greater than 4.3. I presume you meant R >= 4.3

sof202 · 2025-09-05T16:10:12Z

  samples <- names(listSamples)
  for(i in 1:length(listSamples)) {
-    df <-data.frame(listSamples[i])
+    df <-data.frame(listSamples[[i]])


This is good, I encourage indexing like this.

sof202 · 2025-09-05T16:31:50Z

-	multiqc . -f -o ${FASTQCDIR}/multiqc
+	multiqc ${FASTQCDIR}/*_fastqc.zip -f -o ${FASTQCDIR}/multiqc
 	rm -f *.html

-  cd ${ALIGNED_DIR}
-	multiqc . -f -o ${ALIGNED_DIR}/multiqc
+  multiqc ${ALIGNED_DIR}/*.bowtie.log -f -o ${ALIGNED_DIR}/multiqc

-  cd ${TRIM_DIR}/qc
-	multiqc . -f -o ${TRIM_DIR}/multiqc 
+	multiqc ${TRIM_DIR}/qc/*_trimmed_fastqc.zip -f -o ${TRIM_DIR}/multiqc 


Much safer.

sof202 · 2025-09-05T16:34:35Z

-if [[ ! $2 =~ "FRAGSIZE" ]] && [[ ! $2 =~ "PEAKS" ]] && [[ ! $2 =~ "COUNTS" ]] && [[ ! $2 == '' ]];
+if [[ ! $2 =~ "PEAKS" ]] && [[ ! $2 =~ "COUNTS" ]] && [[ ! $2 == 'FILT' ]] && [[ ! $2 == 'CTCHECK' ]] && [[ ! $2 == 'DIFFCOUNTS' ]] && [[ ! $2 == ' ' ]];
 then 
    { echo "Unknown step specified. Please use FRAGSIZE, PEAKS, COUNTS or some combination of this as a single string (i.e. FASTQC,TRIM)" ; exit 1; }            


This is echo has gone stale (COUNTS no longer exists and is replaced by 3 new steps".

If you still want to have combinations of steps you will need to continue using =~ pattern matching instead of == equality.

sof202 · 2025-09-05T16:37:24Z


    ## mark duplicates only
-    picard MarkDuplicates --INPUT ${ALIGNED_DIR}/${sampleName}_sorted_name.bam --OUTPUT ${ALIGNED_DIR}/baseRecalibrate/${sampleName}_dedup.bam --METRICS_FILE ${ALIGNED_DIR}/baseRecalibrate/${sampleName}_metrics.txt --TMP_DIR ${TMPDIR}
+    picard MarkDuplicates --INPUT ${ALIGNED_DIR}/${sampleName}_sorted_name.bam \


Am I correct in assuming that picard creates ${ALIGNED_DIR}/baseRecalibrate? Or is this now created elsewhere (hence the mkdir command being removed on line 37)

marinafloresp and others added 17 commits December 9, 2024 16:32

Add Stage 3 QC

a985bc9

Update README.md

d47d907

Add progress report stage 2 QC

9df4434

Add stage 3 QC scripts

b24e99a

Merge remote-tracking branch 'origin/atac-update' into atac-update

f048a5e

Add stage 3 QC scripts

3f2ccf4

Supress R libraries loading messages

86bb2c4

Rmarkdown lines

13a396c

Rmarkdown lines

6387208

Allow user define peak calling cutoff

6ddf3cd

Update Readme

6db4ba0

Add R libraries

0b7f42c

Minor updates

3a9d610

Update trimming step

80d71c6

Add example file with sequences to trim

e149b9b

Update trimming files

d4e3262

Latest updated ATAC-seq pipeline

329c3b4

marinafloresp added the ATACSeq ATAC-seq data label Sep 5, 2025

sof202 approved these changes Sep 5, 2025

View reviewed changes

Navigation Menu

Update ATAC-seq pipeline - #322

Update ATAC-seq pipeline#322
marinafloresp wants to merge 17 commits into
masterfrom
atac-update

marinafloresp commented Sep 5, 2025 •

edited

Loading

Uh oh!

sof202 left a comment

Uh oh!

sof202 Sep 5, 2025

Uh oh!

sof202 Sep 5, 2025

Uh oh!

sof202 Sep 5, 2025

Uh oh!

sof202 Sep 5, 2025

Uh oh!

sof202 Sep 5, 2025

Uh oh!

sof202 Sep 5, 2025

Uh oh!

sof202 Sep 5, 2025

Uh oh!

sof202 Sep 5, 2025

Uh oh!

sof202 Sep 5, 2025

Uh oh!

sof202 Sep 5, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

		@@ -81,8 +67,6 @@ colors <- brewer.pal(3, "Set2")
		## Load counts in peaks files and match to pheno samples data

Conversation

marinafloresp commented Sep 5, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

Issue ticket number

Type of pull request

Checklist

Uh oh!

sof202 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

marinafloresp commented Sep 5, 2025 •

edited

Loading