Nucleosome phasing

00_Download_and_Preprocessing/ChIP-exo-v6_merge.sh

@@ -49,11 +49,13 @@ java -jar $PICARD MergeSamFiles -O HepG2_IgG_BX_merge_hg38.bam \
 # Use Picard to merge resequenced technical replicates, otherwise rename BAM using cp
 
 # Input
-cp 28326_Input_-_K562_-_-_50Unuclease10min-0cycSonic-newQuenchingbuffer_BI.bam K562_-_BI_rep1_hg38.bam
-
+cp 41859_Input_-_SEM_-_RPMI_-_BI.bam SEM_-_BI_WT_hg38.bam
+cp 41860_Input_-_SEM_-_RPMI_1µM-5PhIAA_BI.bam SEM_-_BI_CTCFKO_hg38.bam
 # CTCF
 cp 33924_CTCF_07-729_K562_-_IMDM_-_BX.bam K562_CTCF_BX_rep1_hg38.bam
 cp 34174_CTCF_07-729_K562_-_IMDM_-_BX.bam K562_CTCF_BX_rep2_hg38.bam
+cp 41780_CTCF_07-729_SEM_-_RPMI_-_BX.bam SEM_CTCF_BX_WT_hg38.bam
+cp 41781_CTCF_07-729_SEM_-_RPMI_CTCF depletion_BX.bam SEM_CTCF_BX_CTCTKO_hg38.bam
 
 # non-crosslinked CTCF 
 java -jar $PICARD MergeSamFiles -O K562_CTCF_nonXLBX_rep1_hg38.bam \
@@ -98,12 +100,10 @@ cp 34158_HNF4A_HPA004712_HepG2_-_-_-_BX.bam HepG2_HNF4A_BX_rep1_hg38.bam
 cp 32659_FOXA2_ab256493_K562_-_-_-_BX.bam K562_FOXA2_BX_rep1_hg38.bam
 cp 38479_FOXA2_ab256493_K562_-_IMDM_-_BX.bam K562_FOXA2_BX_rep2_hg38.bam
 
-# ZKSCAN1 K562
-cp  34048_ZKSCAN1_HPA006672_K562_-_IMDM_-_BX.bam K562_ZKSCAN1_BX_rep1_hg38.bam
-
 cd $WRK/../data
 mv sample-BAM/K562_*.bam $BAMDIR
 mv sample-BAM/HepG2_*.bam $BAMDIR
+mv sample-BAM/SEM_*.bam $BAMDIR
 
 # Index set of BAM files
 for FILE in $BAMDIR/*.bam;

00_Download_and_Preprocessing/X_get_scaling_factors.sbatch

@@ -45,6 +45,7 @@ BAMFILE=`ls $BAMDIR/*hg38.bam | head -n $SLURM_ARRAY_TASK_ID | tail -1`
 BAM=`basename $BAMFILE ".bam"`
 [ -f $BAMFILE.bai ] || samtools index $BAMFILE
 
+
 # Use different normalization method depending on target/assay
 echo "Calculate Total Tag normalization factors w/ blacklist"
 java -jar -Djava.awt.headless=true $SCRIPTMANAGER read-analysis scaling-factor $BAMFILE -f $BLACKLIST --total-tag -o $NDIR/${BAM}_TotalTag
@@ -58,11 +59,15 @@ TARGET=`echo $BAM | cut -d "_" -f 2`
 [ "$STRAIN" == "MNase-ChIP" ] && exit
 [ "$STRAIN" == "CUTRUN" ] && exit
 
+# Calculate NCIS normalization factor and Total tag normalization,  Calculate NCIS normalization factor
+
 CONTROL=$BAMDIR/${STRAIN}_IgG_BX_merge_hg38.bam
 echo "Calculate NCIS normalization factors w/ blacklist and control: ${CONTROL}"
 
 # Calculate NCIS normalization factor
 java -jar -Djava.awt.headless=true $SCRIPTMANAGER read-analysis scaling-factor $BAMFILE -f $BLACKLIST --ncis -c $CONTROL -w 500 -o $NDIR/${BAM}_NCISb
 
 # Calculate fragment size distribution
-java -jar -Djava.awt.headless=true $SCRIPTMANAGER bam-statistics pe-stat -x 500 -d $BAMFILE -o $FDIR/$BAM
+java -jar -Djava.awt.headless=true $SCRIPTMANAGER bam-statistics pe-stat -x 500 -d $BAMFILE -o $FDIR/$BAM
+
+

00_Download_and_Preprocessing/sampleIDs.txt

@@ -31,3 +31,7 @@
 29009
 29010
 34158
+41780
+41781
+41859
+41860

03_Call_JASPAR/2d_SEM_CTCF.sh

@@ -0,0 +1,73 @@
+#!/bin/bash
+
+# data/RefPT-JASPAR-nonK562
+#   |--<TFNAME>_<JASPARID>_Unbound.bed
+#   |--<TFNAME>_<JASPARID>_K562-specific-Unbound.bed
+
+### CHANGE ME
+WRK=/storage/group/bfp2/default/hxc585_HainingChen/2025_Chen_TF-Nuc/03_Call_JASPAR
+#WRK=/ocean/projects/see180003p/owlang/2024-Krebs_Science/03_Call_JASPAR
+#WRK=/scratch/owl5022/2024-Krebs_Science/03_Call_JASPAR
+
+# - java
+
+#set -exo
+module load bedtools
+module load anaconda3
+source activate /storage/group/bfp2/default/owl5022-OliviaLang/conda/bx
+
+# Inputs and outputs
+GENOME=$WRK/../data/hg38_files/hg38.fa
+BLACKLIST=$WRK/../data/hg38_files/ENCFF356LFX_hg38_exclude.bed
+WT_SEM_CTCFBAMFILE=$WRK/../data/BAM/SEM_CTCF_BX_WT_hg38.bam
+CTCFKD_SEM_CTCFBAMFILE=$WRK/../data/BAM/SEM_CTCF_BX_CTCTKO_hg38.bam
+CTCFBEDFOLDER=$WRK/../data/RefPT-JASPAR/
+OUTDIR=$WRK/../data/RefPT-JASPAR-nonK562/
+
+
+# Script shortcuts
+SCRIPTMANAGER=$WRK/../bin/ScriptManager-v0.15.jar
+
+mkdir -p $OUTDIR
+mkdir -p $OUTDIR/1000bp
+mkdir -p $OUTDIR/150bp
+## sort CTCF_MA1929.1.bed by WT_SEM_CTCF, and sort by ratio of CTCFKD_SEM_CTCFBAMFILE to WT_SEM_CTCF
+cat $CTCFBEDFOLDER/CTCF_MA1929.1_SORT-TFnucRatio_GROUP-Quartile*.bed | bedtools sort -i | uniq > $OUTDIR/CTCF_MA1929.1.bed
+java -jar $SCRIPTMANAGER coordinate-manipulation expand-bed -c 100 $OUTDIR/CTCF_MA1929.1.bed -o $OUTDIR/CTCF_MA1929.1_100bp.bed
+java -jar "$SCRIPTMANAGER" read-analysis tag-pileup  $OUTDIR/CTCF_MA1929.1_100bp.bed $WT_SEM_CTCFBAMFILE -1 -s 6  --combined --cpu 4 -M  $OUTDIR/WT_SEM_CTCF_CTCF_MA1929.1_100bp_read1
+java -jar "$SCRIPTMANAGER" read-analysis aggregate-data --sum $OUTDIR/WT_SEM_CTCF_CTCF_MA1929.1_100bp_read1_combined.cdt -o $OUTDIR/WT_SEM_CTCF_CTCF_MA1929.1_100bp_read1_combined.out
+java -jar "$SCRIPTMANAGER" read-analysis tag-pileup  $OUTDIR/CTCF_MA1929.1_100bp.bed $CTCFKD_SEM_CTCFBAMFILE -1 -s 6 --combined --cpu 4 -M  $OUTDIR/CTCTKD_SEM_CTCF_CTCF_MA1929.1_100bp_read1
+java -jar "$SCRIPTMANAGER" read-analysis aggregate-data --sum $OUTDIR/CTCTKD_SEM_CTCF_CTCF_MA1929.1_100bp_read1_combined.cdt -o $OUTDIR/CTCTKD_SEM_CTCF_CTCF_MA1929.1_100bp_read1_combined.out
+
+
+cut -f 2  $OUTDIR/WT_SEM_CTCF_CTCF_MA1929.1_100bp_read1_combined.out | tail -n +2  | cut -f 2 | \
+paste $OUTDIR/CTCF_MA1929.1.bed - > $OUTDIR/CTCF_MA1929.1_WT-SEM-CTCF.tsv
+cut -f 2  $OUTDIR/CTCTKD_SEM_CTCF_CTCF_MA1929.1_100bp_read1_combined.out | tail -n +2  | cut -f 2 | \
+paste $OUTDIR/CTCF_MA1929.1_WT-SEM-CTCF.tsv - > $OUTDIR/CTCF_MA1929.1_WT-SEM-CTCF_CTCTKD_SEM_CTCF.tsv
+
+rm  $OUTDIR/WT_SEM_CTCF_CTCF_MA1929.1_100bp_read1_combined.cdt $OUTDIR/CTCTKD_SEM_CTCF_CTCF_MA1929.1_100bp_read1_combined.cdt
+rm   $OUTDIR/CTCTKD_SEM_CTCF_CTCF_MA1929.1_100bp_read1_combined.out $OUTDIR/WT_SEM_CTCF_CTCF_MA1929.1_100bp_read1_combined.out
+rm $OUTDIR/CTCF_MA1929.1_100bp.bed
+rm  $OUTDIR/CTCF_MA1929.1_WT-SEM-CTCF.tsv
+
+
+bedtools sort -i $OUTDIR/CTCF_MA1929.1_WT-SEM-CTCF_CTCTKD_SEM_CTCF.tsv  | awk -v DIR="$OUTDIR" 'BEGIN{OFS="\t";FS="\t"}{
+    if ($7 >= 1) {
+        system("echo \"" $0 "\" >> " DIR "/CTCF_MA1929.1_SORT-WT-SEM-CTCF.tsv")
+    }
+}'  
+
+awk  'BEGIN{FS="\t";OFS="\t"}{print $1,$2,$3,$4,$5,$6,$7,$8,$8/$7}' $OUTDIR/CTCF_MA1929.1_SORT-WT-SEM-CTCF.tsv | sort -k9,9nr > $OUTDIR/CTCF_MA1929.1_SORT-CTCFKD-to-WT-SEM-CTCF.tsv
+
+rm $OUTDIR/CTCF_MA1929.1_SORT-WT-SEM-CTCF.tsv
+
+wc -l $OUTDIR/CTCF_MA1929.1_SORT-CTCFKD-to-WT-SEM-CTCF.tsv
+#  16671 /CTCF_MA1929.1_SORT-CTCFKD-to-WT-SEM-CTCF.tsv
+
+sort -k7,7nr $OUTDIR/CTCF_MA1929.1_SORT-CTCFKD-to-WT-SEM-CTCF.tsv  | head -n 4168 | cut -f 1-6 - > $OUTDIR/CTCF_MA1929.1_1.bed
+
+for file in  $OUTDIR/CTCF_MA1929.1_1.bed  ; do
+	filename=`basename $file ".bed"`
+	java -jar $SCRIPTMANAGER coordinate-manipulation expand-bed -c 1000 $file -o $OUTDIR/1000bp/${filename}_1000bp.bed
+	java -jar $SCRIPTMANAGER coordinate-manipulation expand-bed -c 150 $file -o $OUTDIR/150bp/${filename}_150bp.bed
+done

03_Call_JASPAR/README.md

@@ -119,5 +119,18 @@ Along with intermediate files, two final motif RefPT files are generated for eac
 ../data/RefPT-JASPAR-nonK562/MA1929_1_mm10_intersected_MPE-seq10min_164bp_category*.bed
 ../data/RefPT-JASPAR-nonK562/1000bp/MA1929_1_mm10_intersected_MPE-seq10min_164bp_category*_1000bp.bed
 ```
-
+### 2d_SEM_CTCF.sh
+Get CTCF motif binding sites from SEM cell treate with DMSO or 1µM 5PhIAA 
+```
+sh 2d_SEM_CTCF.sh
+```
+Along with intermediate files, three final motif RefPT files are generated for CTCF 4 group depleted or undepleted:
+```
+../data/RefPT-JASPAR-nonK562/CTCF_MA1929.1_*_KO-nondepleted.bed
+../data/RefPT-JASPAR-nonK562/1000bp/CTCF_MA1929.1_*_KO-nondepleted_1000bp.bed
+../data/RefPT-JASPAR-nonK562/150bp/CTCF_MA1929.1_*_KO-nondepleted_150bp.bed
+../data/RefPT-JASPAR-nonK562/CTCF_MA1929.1_*_KO-depleted.bed
+../data/RefPT-JASPAR-nonK562/1000bp/CTCF_MA1929.1_*_KO-depleted_1000bp.bed
+../data/RefPT-JASPAR-nonK562/150bp/CTCF_MA1929.1_*_KO-depleted_150bp.bed
+```
 ```

README.md

@@ -1,21 +1,19 @@
 # Genome-wide rotational and translational setting of transcription factors with nucleosomes
 
-### Jordan E. Krebs<sup>&#x2020; 1,2</sup>, Haining Chen<sup>&#x2020; 2</sup>, Olivia W. Lang<sup>2</sup>, William K. M. Lai<sup>2</sup>, B. Franklin Pugh<sup>2</sup>
+### Haining Chen<sup>&#x2020; 1</sup>, Jordan E. Krebs<sup>&#x2020; 1</sup>,  Olivia W. Lang<sup>1</sup>, Judith H. Lang<sup>1</sup>, Chunliang L. Lai<sup>1</sup>, William K. M. Lai<sup>1</sup>, B. Franklin Pugh<sup>1</sup>
 
 &#x2020; Co-first author
 
-<sup>1</sup>MD/PhD Medical Scientist Training Program, Penn State College of Medicine, Hershey, PA, USA.
-
-<sup>2</sup>Department of Molecular Biology and Genetics, Cornell University, Ithaca, New York, 14853, USA
-
+<sup>1</sup>Department of Molecular Biology and Genetics, Cornell University, Ithaca, New York, 14853, USA
+<sup>1</sup>Department of Tumor Cell Biology, St. Jude Children’s Research Hospital, Memphis, Tennessee, 38105, USA
 ### Correspondence:fp265@cornell.edu
 
 ### PMID : [XXXXXXXX](https://pubmed.ncbi.nlm.nih.gov/XXXXXXXX/)
 ### GEO ID : [GSE266547](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE266547)
+### GEO ID : [GSE267711](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE267711)
 
 ## Abstract
-How transcription factors (TFs) access DNA that is packed into chromatin has been challenging to decipher on a genomic scale due to inherent resolution limits of assays. Accessibility differs depending on whether a TF binding site (TFBS) is translationally and/or rotationally phased on or next to a nucleosome. Rotational phasing occurs where a DNA sequence in the DNA helix has a predominant rotational direction on the nucleosome surface, thereby consistently facing outward (accessible) or inward (inaccessible). While rotational phasing is DNA-encoded by dinucleotide periodicities at yeast TFBSs, such encoding has not been found in humans. Here we develop a genome-wide Benzonase nuclease-based assay to measure translational and rotational phasing, and a second assay to measure such phasing on the same DNA molecule to which a TF is bound. The latter uniquely allows phasing to be measured with transient TFs. We show that many types of human TFBSs have distinct translational and/or rotational phasing depending on whether they are bound by specific TFs. For example, unbound CTCF sites are nucleosomal with local but not global translational phasing. They possess DNA-encoded rotational phasing, with the predominant phase being in an accessible orientation. CTCF-bound sites have adjacent nucleosomes possessing a global translational phase and distinct DNA-encoded rotational phases. Similar themes recur for other TFs including NFIA and pioneer factor FoxA, but their activities are highly transient. Our assays further allow subnucleosomal structures to be examined along with their relationship to a transiting RNA polymerase. Together, these findings reveal an intimate relationship between nucleosome phasing and transcription factors.
-
+How transcription factors (TFs) access DNA and interface within nucleosomes is key to understanding gene regulation. TFs recognize their DNA binding site (TFBS) either when it is nucleosome-free or its helical site is rotationally accessible on a nucleosome surface. These two access parameters are manifested through translational and rotational phasing of nucleosomes. Measuring phasing on a genomic scale has been limited by the resolution of current assays. Here we develop a nuclease-based assay (Benzonase-seq) to probe translational and rotational phasing across the human genome, and update ChIP-exo (to v6) to measure this phasing on the same DNA molecule to which a TF is bound. The latter uniquely allows phasing to be measured with TFs that bind DNA transiently. We find that human TFBSs show distinct phasing patterns depending on TF occupancy. For example, unbound CTCF sites are encased within DNA-encoded rotationally phased nucleosomes. But there is no consensus translational phasing until CTCF binds. Similar relationships exist with other TFs like ZKSCAN1, and for transiently binding TFs like FoxA and NFIA. Our assays also provide improved maps of nucleosome positions particularly at promoters and resolves subnucleosomal structures. These findings uncover new principles of how TFs engage nucleosomes to shape the genomic chromatin landscape.
 
 ## Directions
 To recreate the figures for this manuscript, please execute the scripts in each directory in numerical order. Each directory's README includes more specific details on execution. To be more explicit, run the scripts in each directory in the following order: `00_Download_and_Preprocessing`, `01_Run_GenoPipe`, `02_Call_Nucleosomes`, `03_Call_JASPAR`, `04_Call_Motifs`, `X_Bulk_Processing`, and then finally `Z_Figures`.

X_Bulk_Processing/8d_CTCF_Q4_CTCF-KO.sbatch

@@ -0,0 +1,74 @@
+#!/bin/bash
+
+# CTCF KD vs CTCF control BI
+
+### CHANGE ME
+WRK=/storage/group/bfp2/default/hxc585_HainingChen/2025_Chen_TF-Nuc/X_Bulk_Processing
+
+###
+
+# Dependencies
+# - bedtools
+# - java
+
+set -exo
+module load bedtools
+module load anaconda3
+source activate /storage/group/bfp2/default/owl5022-OliviaLang/conda/bx
+
+# Inputs and outputs
+GENOME=$WRK/../data/hg38_files/hg38.fa
+Genome=$WRK/../data/hg38_files/hg38.info.txt
+BAMDIR=$WRK/../data/BAM
+BICTCFKD=$WRK/../data/BAM/SEM_-_BI_CTCFKO_hg38.bam
+BICTCFControl=$WRK/../data/BAM/SEM_-_BI_WT_hg38.bam
+WT_SEM_CTCFBAMFILE=$WRK/../data/BAM/SEM_CTCF_BX_WT_hg38.bam
+CTCFKD_SEM_CTCFBAMFILE=$WRK/../data/BAM/SEM_CTCF_BX_CTCTKO_hg38.bam
+BLACKLIST=$WRK/../data/hg38_files/ENCFF356LFX_hg38_exclude.bed
+BEDDIR=$WRK/../data/RefPT-JASPAR-nonK562/
+DIR=$WRK/Library/CTCF_MA1929.1
+
+# cd to data file
+cd $WRK/
+# Script shortcuts
+
+SCRIPTMANAGER=$WRK/../bin/ScriptManager-v0.15.jar
+COMPOSITE=$WRK/../bin/sum_Col_CDT.pl
+
+## Create output directories
+[ -d $DIR ] || mkdir $DIR
+[[ -d $DIR/SCORES ]] || mkdir $DIR/SCORES
+
+## pile up
+
+for file in $BICTCFKD  $BICTCFControl ; do
+    BAM=`basename $file ".bam"`
+    java -jar "$SCRIPTMANAGER"  bam-statistics se-stat $file -o $DIR/SCORES/${BAM}.out
+done
+
+for file in $WT_SEM_CTCFBAMFILE  $CTCFKD_SEM_CTCFBAMFILE    ; do
+    BAM=`basename $file ".bam"`
+    NFFILE=$BAMDIR/NormalizationFactors/$BAM\_TotalTag_ScalingFactors.out
+    FACTOR=`grep 'Scaling factor' $NFFILE | awk -F" " '{print $3}'`
+    java -jar $SCRIPTMANAGER read-analysis tag-pileup $BEDDIR/1000bp/CTCF_MA1929.1_1_1000bp.bed $file  --cpu 4 -1 -o $DIR/${BAM}_CTCF_MA1929.1_1_1000bp_read1.out
+    java -jar $SCRIPTMANAGER read-analysis scale-matrix $DIR/${BAM}_CTCF_MA1929.1_1_1000bp_read1.out  -s $FACTOR -r 1 -l 1 -o $DIR/${BAM}_CTCF_MA1929.1_1_1000bp_read1_Normalized.out
+    rm $DIR/${BAM}_CTCF_MA1929.1_1_1000bp_read1.out
+    java -jar $SCRIPTMANAGER read-analysis tag-pileup $BEDDIR/150bp/CTCF_MA1929.1_1_150bp.bed $file  --cpu 4 -s 6 --combined  -M $DIR/${BAM}_CTCF_MA1929.1_1_150bp_read1
+    java -jar $SCRIPTMANAGER read-analysis scale-matrix $DIR/${BAM}_CTCF_MA1929.1_1_150bp_read1_combined.cdt  -s $FACTOR -o $DIR/${BAM}_CTCF_MA1929.1_1_150bp_read1_combined_Normalized.out
+    java -jar "$SCRIPTMANAGER" read-analysis aggregate-data --sum $DIR/${BAM}_CTCF_MA1929.1_1_150bp_read1_combined_Normalized.out -o $DIR/SCORES/${BAM}_CTCF_MA1929.1_1_150bp_read1_combined_Normalized.out
+    rm $DIR/${BAM}_CTCF_MA1929.1_1_150bp_read1_combined.cdt 
+done
+
+for file in $BICTCFKD  $BICTCFControl ; do
+    BAM=`basename $file ".bam"`
+    NFFILE=$BAMDIR/NormalizationFactors/$BAM\_TotalTag_ScalingFactors.out
+    FACTOR=`grep 'Scaling factor' $NFFILE | awk -F" " '{print $3}'`
+    java -jar $SCRIPTMANAGER read-analysis tag-pileup $BEDDIR/1000bp/CTCF_MA1929.1_1_1000bp.bed $file  --cpu 4 -m --combined -o $DIR/${BAM}_CTCF_MA1929.1_1_1000bp_midpoint.out 
+    java -jar $SCRIPTMANAGER read-analysis scale-matrix $DIR/${BAM}_CTCF_MA1929.1_1_1000bp_midpoint.out  -s $FACTOR -r 1 -l 1 -o $DIR/${BAM}_CTCF_MA1929.1_1_1000bp_midpoint_Normalized.cdt
+    rm $DIR/${BAM}_CTCF_MA1929.1_1_1000bp_midpoint.out
+    java -jar $SCRIPTMANAGER read-analysis tag-pileup $BEDDIR/150bp/CTCF_MA1929.1_1_150bp.bed $file  --cpu 4 -m --combined  -M $DIR/${BAM}_CTCF_MA1929.1_1_150bp_midpoint
+    java -jar "$SCRIPTMANAGER" read-analysis aggregate-data --sum $DIR/${BAM}_CTCF_MA1929.1_1_150bp_midpoint_combined.cdt -o $DIR/SCORES/${BAM}_CTCF_MA1929.1_1_150bp_midpoint_combined.out
+    rm $DIR/${BAM}_CTCF_MA1929.1_1_150bp_midpoint_combined.cdt 
+done
+
+