# ATAC-Seq

## mapping

for input in `seq 1 16`

do

echo "-----starting ${input} mapping-----" \

gunzip ./fastq/sample${input}_S${input}_R1_001.fastq.gz
gunzip ./fastq/sample${input}_S${input}_R2_001.fastq.gz

bowtie2 -x /home/t-t/bowtie2_index_mm10_ensmbl_chr/bowtie2_index_mm10_ensmbl_chr -1 ./fastq/sample${input}_S${input}_R1_001.fastq -2 ./fastq/sample${input}_S${input}_R2_001.fastq -S ./mapping/sample${input}.sam -q -N 1 -p 10 2> sample${input}.log

samtools view -@ 10 -bS ./mapping/sample${input}.sam -o ./bam/sample${input}.bam

samtools sort -@ 10 ./bam/sample${input}.bam -o ./bam/sample${input}_sorted.bam

samtools index ./bam/sample${input}_sorted.bam

done

## duplicate removal

for input in `seq 1 16`

do
picard MarkDuplicates I=sample${input}_sorted.bam O=sample${input}_sorted_rmdup.bam M=sample${input}_sorted_rmdup.bam.log.txt REMOVE_DUPLICATES=true

samtools index sample${input}_sorted_rmdup.bam
done

## フラグメント長（インサートサイズ）の分布の可視化

for input in `seq 1 16`

do

picard CollectInsertSizeMetrics INPUT=bam/sample${input}_sorted_rmdup.bam OUTPUT=picard/sample${input}_insert_size_metrics.txt HISTOGRAM_FILE=picard/sample${input｝_hist.pdf MINIMUM_PCT=0

done

## deeptools

for input in `seq 1 16`

do

bamCoverage -bs 5 -p max --normalizeUsing RPKM --effectiveGenomeSize 2913022398 -b ./bam/sample${input}_sorted_rmdup.bam -o ./bigwig/sample${input}.bw

done

## peakcall

for input in `seq 1 16`

do

macs2 callpeak -t bam/sample${input}_sorted_rmdup.bam -f BAM --nomodel --shift -50 --extsize 100 -g mm -n sample${input} -B --outdir macs2

done

#!/bin/bash

for input in `seq 1 16`

do

cat 230714_ATAC/fastq/sample${input}_S${input}_R1_001.fastq 230930_ATAC_re/fastq/sample${input}_S${input}_R1_001.fastq > 230930_ATAC_re/fastq_merge/sample${input}_R1.fastq

cat 230714_ATAC/fastq/sample${input}_S${input}_R2_001.fastq 230930_ATAC_re/fastq/sample${input}_S${input}_R2_001.fastq > 230930_ATAC_re/fastq_merge/sample${input}_R2.fastq

done

ATAC-Seq参考
https://informatics.fas.harvard.edu/atac-seq-guidelines.html
http://barc.wi.mit.edu/education/hot_topics/ATACseq_2023/ATACseq2023_4slidesPerPage.pdf

###

bowtie2 -x /home/t-t/bowtie2_index_mm10_ensmbl_chr/bowtie2_index_mm10_ensmbl_chr -1 ./fastq_merge/sample${input}_R1.fastq -2 ./fastq_merge/sample${input}_R2.fastq -S ./mapping/sample${input}.sam -q -N 1 -p 10 2> sample${input}.log

samtools view -@ 10 -bS ./mapping/sample${input}.sam -o ./bam/sample${input}.bam

samtools sort -@ 10 ./bam/sample${input}.bam -o ./bam/sample${input}_sorted.bam

picard MarkDuplicates I=./bam/sample${input}_sorted.bam O=./bam/sample${input}_sorted_rmdup.bam M=./bam/sample${input}_sorted_rmdup.bam.log.txt REMOVE_DUPLICATES=true

samtools index ./bam/sample${input}_sorted_rmdup.bam

bamCoverage -bs 5 -p max --normalizeUsing RPKM --effectiveGenomeSize 2913022398 -b ./bam/sample${input}_sorted_rmdup.bam -o ./bigwig/sample${input}.bw

macs2 callpeak -t bam/sample${input}_sorted_rmdup.bam -f BAM --nomodel --shift -50 --extsize 100 -g mm -n sample${input} -B --outdir macs2

### remove reads mapped to mitochondria

samtools view -@ 8 -f 0x2 -b ./bam/sample${input}_sorted_rmdup.bam chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chrX chrY > ./bam/chrM/sample${input}_sorted_rmdup_chrM.bam

samtools index ./bam/chrM/sample${input}_sorted_rmdup_chrM.bam

bamCoverage -bs 5 -p max --normalizeUsing CPM --effectiveGenomeSize 2913022398 -b ./bam/chrM/sample${input}_sorted_rmdup_chrM.bam -o ./bigwig/chrM/sample${input}.bw

macs2 callpeak -t bam/chrM/sample${input}_sorted_rmdup_chrM.bam -f BAM --nomodel --shift -50 --extsize 100 -g mm -n sample${input} -B --outdir macs2/chrM/

N=2サンプルのオーバーラップした領域を抽出 (bedtools intersect)

bedtools intersect -a macs/sample1_peaks.narrowPeak -b macs/sample5_peaks.narrowPeak -sorted > macs/overlap/peaks1_5.bed
bedtools intersect -a macs/sample9_peaks.narrowPeak -b macs/sample13_peaks.narrowPeak -sorted > macs/overlap/peaks9_13.bed
bedtools intersect -a macs/sample2_peaks.narrowPeak -b macs/sample6_peaks.narrowPeak -sorted > macs/overlap/peaks2_6.bed
bedtools intersect -a macs/sample10_peaks.narrowPeak -b macs/sample14_peaks.narrowPeak -sorted > macs/overlap/peaks10_14.bed
bedtools intersect -a macs/sample3_peaks.narrowPeak -b macs/sample7_peaks.narrowPeak -sorted > macs/overlap/peaks3_7.bed
bedtools intersect -a macs/sample11_peaks.narrowPeak -b macs/sample15_peaks.narrowPeak -sorted > macs/overlap/peaks11_15.bed
bedtools intersect -a macs/sample4_peaks.narrowPeak -b macs/sample8_peaks.narrowPeak -sorted > macs/overlap/peaks4_8.bed
bedtools intersect -a macs/sample12_peaks.narrowPeak -b macs/sample16_peaks.narrowPeak -sorted > macs/overlap/peaks12_16.bed

オーバーラップしない領域を抽出　(bedtools intersect -v)
WT特異的な領域

bedtools intersect -v -a macs/overlap/peaks1_5.bed -b macs/overlap/peaks9_13.bed -sorted > macs/specific/DP_WT.bed
bedtools intersect -v -a macs/overlap/peaks2_6.bed -b macs/overlap/peaks10_14.bed -sorted > macs/specific/pre1_WT.bed
bedtools intersect -v -a macs/overlap/peaks3_7.bed -b macs/overlap/peaks11_15.bed -sorted > macs/specific/pre2_WT.bed
bedtools intersect -v -a macs/overlap/peaks4_8.bed -b macs/overlap/peaks12_16.bed -sorted > macs/specific/mature_WT.bed

KO特異的な領域

bedtools intersect -v -b macs/overlap/peaks1_5.bed -a macs/overlap/peaks9_13.bed -sorted > macs/specific/DP_KO.bed
bedtools intersect -v -b macs/overlap/peaks2_6.bed -a macs/overlap/peaks10_14.bed -sorted > macs/specific/pre1_KO.bed
bedtools intersect -v -b macs/overlap/peaks3_7.bed -a macs/overlap/peaks11_15.bed -sorted > macs/specific/pre2_KO.bed
bedtools intersect -v -b macs/overlap/peaks4_8.bed -a macs/overlap/peaks12_16.bed -sorted > macs/specific/mature_KO.bed

DP -> Pre1 -> pre2 -> matureで特異的な遺伝子座

# DPの中でpre1にない（DP特異的）
bedtools intersect -v -a macs/overlap/peaks1_5.bed -b macs/overlap/peaks2_6.bed -sorted > macs/specific/DP_specific.bed
# pre1の中でDPにない（pre1特異的）
bedtools intersect -v -b macs/overlap/peaks1_5.bed -a macs/overlap/peaks2_6.bed -sorted > macs/specific/pre1a_specific.bed

# pre1の中でpre2にない（pre1特異的）
bedtools intersect -v -a macs/overlap/peaks2_6.bed -b macs/overlap/peaks3_7.bed -sorted > macs/specific/pre1b_specific.bed

# pre2の中でpre1にない（pre2特異的）
bedtools intersect -v -b macs/overlap/peaks2_6.bed -a macs/overlap/peaks3_7.bed -sorted > macs/specific/pre2a_specific.bed

# pre2の中でmatureにない（pre2特異的）
bedtools intersect -v -a macs/overlap/peaks3_7.bed -b macs/overlap/peaks4_8.bed -sorted > macs/specific/pre2b_specific.bed

# matureの中でpre2にない（mature特異的）
bedtools intersect -v -b macs/overlap/peaks3_7.bed -a macs/overlap/peaks4_8.bed -sorted > macs/specific/mature_specific.bed

N=2サンプルピークをマージする (bedtools merge)
100bp以内のピーク同士はmergeする (-d 100オプション)

cat macs/sample1_peaks.narrowPeak macs/sample5_peaks.narrowPeak | sort -k1,1 -k2,2n | bedtools merge -d 100 > macs/merge/peak1_5_merge.bed
cat macs/sample9_peaks.narrowPeak macs/sample13_peaks.narrowPeak | sort -k1,1 -k2,2n | bedtools merge -d 100 > macs/merge/peak9_13_merge.bed

cat macs/sample2_peaks.narrowPeak macs/sample6_peaks.narrowPeak | sort -k1,1 -k2,2n | bedtools merge -d 100 > macs/merge/peak2_6_merge.bed

cat macs/sample10_peaks.narrowPeak macs/sample14_peaks.narrowPeak | sort -k1,1 -k2,2n | bedtools merge -d 100 > macs/merge/peak10_14_merge.bed

cat macs/sample3_peaks.narrowPeak macs/sample7_peaks.narrowPeak | sort -k1,1 -k2,2n | bedtools merge -d 100 > macs/merge/peak3_7_merge.bed

cat macs/sample11_peaks.narrowPeak macs/sample15_peaks.narrowPeak | sort -k1,1 -k2,2n | bedtools merge -d 100 > macs/merge/peak11_15_merge.bed

cat macs/sample4_peaks.narrowPeak macs/sample8_peaks.narrowPeak | sort -k1,1 -k2,2n | bedtools merge -d 100 > macs/merge/peak4_8_merge.bed

cat macs/sample12_peaks.narrowPeak macs/sample16_peaks.narrowPeak | sort -k1,1 -k2,2n | bedtools merge -d 100 > macs/merge/peak12_16_merge.bed

## ChipPeakAnnoのインストール

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("ChIPpeakAnno")

library(ChIPpeakAnno)

データの格納

# 全ピークコール領域
DP_WT_all <- toGRanges("overlap/peaks1_5.bed", format="narrowPeak", header = F)
pre1_WT_all <- toGRanges("overlap/peaks2_6.bed", format="narrowPeak", header = F)
pre2_WT_all <- toGRanges("overlap/peaks3_7.bed", format="narrowPeak", header = F)
mature_WT_all <- toGRanges("overlap/peaks4_8.bed", format="narrowPeak", header = F)
DP_KO_all <- toGRanges("overlap/peaks9_13.bed", format="narrowPeak", header = F)
pre1_KO_all <- toGRanges("overlap/peaks10_14.bed", format="narrowPeak", header = F)
pre2_KO_all <- toGRanges("overlap/peaks11_15.bed", format="narrowPeak", header = F)
mature_KO_all <- toGRanges("overlap/peaks12_16.bed", format="narrowPeak", header = F)

# WT特異的
DP_WT <- toGRanges("specific/DP_WT.bed", format="narrowPeak", header = F)
pre1_WT <- toGRanges("specific/pre1_WT.bed", format="narrowPeak", header = F)
pre2_WT <- toGRanges("specific/pre2_WT.bed", format="narrowPeak", header = F)
mature_WT <- toGRanges("specific/mature_WT.bed", format="narrowPeak", header = F)

# KO特異的
DP_KO <- toGRanges("specific/DP_KO.bed", format="narrowPeak", header = F)
pre1_KO <- toGRanges("specific/pre1_KO.bed", format="narrowPeak", header = F)
pre2_KO <- toGRanges("specific/pre2_KO.bed", format="narrowPeak", header = F)
mature_KO <- toGRanges("specific/mature_KO.bed", format="narrowPeak", header = F)

# developmental stage特異的（WTのみ）
DP_specific <- toGRanges("specific/DP_specific.bed", format="narrowPeak", header = F)
pre1a_specific <- toGRanges("specific/pre1a_specific.bed", format="narrowPeak", header = F)
pre1b_specific <- toGRanges("specific/pre1b_specific.bed", format="narrowPeak", header = F)
pre2a_specific <- toGRanges("specific/pre2a_specific.bed", format="narrowPeak", header = F)
pre2b_specific <- toGRanges("specific/pre2b_specific.bed", format="narrowPeak", header = F)
mature_specific <- toGRanges("specific/mature_specific.bed", format="narrowPeak", header = F)

ピークに最も近い転写開始点の遺伝子へ割り当てる
マウスmm10遺伝子モデルのパッケージのダウンロード

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("TxDb.Mmusculus.UCSC.mm10.ensGene")

マウスmm10の遺伝子モデルのダウンロードとロード

if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("TxDb.Mmusculus.UCSC.mm10.ensGene")
library(TxDb.Mmusculus.UCSC.mm10.ensGene)
annoData <- toGRanges(TxDb.Mmusculus.UCSC.mm10.ensGene)

seqlevelsStyle(DP_WT_all) <- seqlevelsStyle(annoData)
seqlevelsStyle(pre1_WT_all) <- seqlevelsStyle(annoData)
seqlevelsStyle(pre2_WT_all) <- seqlevelsStyle(annoData)
seqlevelsStyle(mature_WT_all) <- seqlevelsStyle(annoData)
seqlevelsStyle(DP_KO_all) <- seqlevelsStyle(annoData)
seqlevelsStyle(pre1_KO_all) <- seqlevelsStyle(annoData)
seqlevelsStyle(pre2_KO_all) <- seqlevelsStyle(annoData)
seqlevelsStyle(mature_KO_all) <- seqlevelsStyle(annoData)
seqlevelsStyle(DP_WT) <- seqlevelsStyle(annoData)
seqlevelsStyle(pre1_WT) <- seqlevelsStyle(annoData)
seqlevelsStyle(pre2_WT) <- seqlevelsStyle(annoData)
seqlevelsStyle(mature_WT) <- seqlevelsStyle(annoData)
seqlevelsStyle(DP_KO) <- seqlevelsStyle(annoData)
seqlevelsStyle(pre1_KO) <- seqlevelsStyle(annoData)
seqlevelsStyle(pre2_KO) <- seqlevelsStyle(annoData)
seqlevelsStyle(mature_KO) <- seqlevelsStyle(annoData)
seqlevelsStyle(DP_specific) <- seqlevelsStyle(annoData)
seqlevelsStyle(pre1a_specific) <- seqlevelsStyle(annoData)
seqlevelsStyle(pre1b_specific) <- seqlevelsStyle(annoData)
seqlevelsStyle(pre2a_specific) <- seqlevelsStyle(annoData)
seqlevelsStyle(pre2b_specific) <- seqlevelsStyle(annoData)
seqlevelsStyle(mature_specific) <- seqlevelsStyle(annoData)

DP_WT_all_anno <- annotatePeakInBatch(DP_WT_all,AnnotationData = annoData)
pre1_WT_all_anno <- annotatePeakInBatch(pre1_WT_all,AnnotationData = annoData)
pre2_WT_all_anno <- annotatePeakInBatch(pre2_WT_all,AnnotationData = annoData)
mature_WT_all_anno <- annotatePeakInBatch(mature_WT_all,AnnotationData = annoData)
DP_KO_all_anno <- annotatePeakInBatch(DP_KO_all,AnnotationData = annoData)
pre1_KO_all_anno <- annotatePeakInBatch(pre1_KO_all,AnnotationData = annoData)
pre2_KO_all_anno <- annotatePeakInBatch(pre2_KO_all,AnnotationData = annoData)
mature_KO_all_anno <- annotatePeakInBatch(mature_KO_all,AnnotationData = annoData)
DP_WT_anno <- annotatePeakInBatch(DP_WT,AnnotationData = annoData)
pre1_WT_anno <- annotatePeakInBatch(pre1_WT,AnnotationData = annoData)
pre2_WT_anno <- annotatePeakInBatch(pre2_WT,AnnotationData = annoData)
mature_WT_anno <- annotatePeakInBatch(mature_WT,AnnotationData = annoData)
DP_KO_anno <- annotatePeakInBatch(DP_KO,AnnotationData = annoData)
pre1_KO_anno <- annotatePeakInBatch(pre1_KO,AnnotationData = annoData)
pre2_KO_anno <- annotatePeakInBatch(pre2_KO,AnnotationData = annoData)
mature_KO_anno <- annotatePeakInBatch(mature_KO,AnnotationData = annoData)

DP_sepcific_anno <- annotatePeakInBatch(DP_specific,AnnotationData = annoData)
pre1a_sepcific_anno <- annotatePeakInBatch(pre1a_specific,AnnotationData = annoData)
pre1b_sepcific_anno <- annotatePeakInBatch(pre1b_specific,AnnotationData = annoData)
pre2a_sepcific_anno <- annotatePeakInBatch(pre2a_specific,AnnotationData = annoData)
pre2b_sepcific_anno <- annotatePeakInBatch(pre2b_specific,AnnotationData = annoData)
mature_sepcific_anno <- annotatePeakInBatch(mature_specific,AnnotationData = annoData)

png("DP_WT_all.png", width = 2000, height =1600, res = 400)
pie1(table(DP_WT_all_anno$insideFeature), main = "DP_WT_all")
dev.off()
png("pre1_WT_all.png", width = 2000, height =1600, res = 400)
pie1(table(pre1_WT_all_anno$insideFeature), main = "pre1_WT_all")
dev.off()
png("pre2_WT_all.png", width = 2000, height =1600, res = 400)
pie1(table(pre2_WT_all_anno$insideFeature), main = "pre2_WT_all")
dev.off()
png("mature_WT_all.png", width = 2000, height =1600, res = 400)
pie1(table(mature_WT_all_anno$insideFeature), main = "mature_WT_all")
dev.off()
png("DP_KO_all.png", width = 2000, height =1600, res = 400)
pie1(table(DP_KO_all_anno$insideFeature), main = "DP_KO_all")
dev.off()
png("pre1_KO_all.png", width = 2000, height =1600, res = 400)
pie1(table(pre1_KO_all_anno$insideFeature), main = "pre1_KO_all")
dev.off()
png("pre2_KO_all.png", width = 2000, height =1600, res = 400)
pie1(table(pre2_KO_all_anno$insideFeature), main = "pre2_KO_all")
dev.off()
png("mature_KO_all.png", width = 2000, height =1600, res = 400)
pie1(table(mature_KO_all_anno$insideFeature), main = "mature_KO_all")
dev.off()
png("DP_WT.png", width = 2000, height =1600, res = 400)
pie1(table(DP_WT_anno$insideFeature), main = "DP_WT")
dev.off()
png("pre1_WT.png", width = 2000, height =1600, res = 400)
pie1(table(pre1_WT_anno$insideFeature), main = "pre1_WT")
dev.off()
png("pre2_WT.png", width = 2000, height =1600, res = 400)
pie1(table(pre2_WT_anno$insideFeature), main = "pre2_WT")
dev.off()
png("mature_WT.png", width = 2000, height =1600, res = 400)
pie1(table(mature_WT_anno$insideFeature), main = "mature_WT")
dev.off()
png("DP_KO.png", width = 2000, height =1600, res = 400)
pie1(table(DP_KO_anno$insideFeature), main = "DP_KO")
dev.off()
png("pre1_KO.png", width = 2000, height =1600, res = 400)
pie1(table(pre1_KO_anno$insideFeature), main = "pre1_KO")
dev.off()
png("pre2_KO.png", width = 2000, height =1600, res = 400)
pie1(table(pre2_KO_anno$insideFeature), main = "pre2_KO")
dev.off()
png("mature_KO.png", width = 2000, height =1600, res = 400)
pie1(table(mature_KO_anno$insideFeature), main = "mature_KO")
dev.off()


png("DP_sepcific.png", width = 2000, height =1600, res = 400)
pie1(table(DP_sepcific_anno$insideFeature), main = "DP_specific")
dev.off()
png("pre1a_specific.png", width = 2000, height =1600, res = 400)
pie1(table(pre1a_sepcific_anno$insideFeature), main = "pre1a_specific")
dev.off()
png("pre1b_specific.png", width = 2000, height =1600, res = 400)
pie1(table(pre1b_sepcific_anno$insideFeature), main = "pre1b_specific")
dev.off()
png("pre2a_specific.png", width = 2000, height =1600, res = 400)
pie1(table(pre2a_sepcific_anno$insideFeature), main = "pre2a_specific")
dev.off()
png("pre2b_specific.png", width = 2000, height =1600, res = 400)
pie1(table(pre2b_sepcific_anno$insideFeature), main = "pre2b_specific")
dev.off()
png("mature_specific.png", width = 2000, height =1600, res = 400)
pie1(table(mature_sepcific_anno$insideFeature), main = "mature_specific")
dev.off()

遺伝子IDに対応する遺伝子名を追加

BiocManager::install("EnsDb.Mmusculus.v79")
library(EnsDb.Mmusculus.v79)

DP_WT_all_anno$feature[is.na(DP_WT_all_anno$feature)] <- "."
DP_WT_all_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = DP_WT_all_anno$feature, column = "GENENAME", keytype = "GENEID")
pre1_WT_all_anno$feature[is.na(pre1_WT_all_anno$feature)] <- "."
pre1_WT_all_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = pre1_WT_all_anno$feature, column = "GENENAME", keytype = "GENEID")
pre2_WT_all_anno$feature[is.na(pre2_WT_all_anno$feature)] <- "."
pre2_WT_all_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = pre2_WT_all_anno$feature, column = "GENENAME", keytype = "GENEID")
mature_WT_all_anno$feature[is.na(mature_WT_all_anno$feature)] <- "."
mature_WT_all_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = mature_WT_all_anno$feature, column = "GENENAME", keytype = "GENEID")
DP_KO_all_anno$feature[is.na(DP_KO_all_anno$feature)] <- "."
DP_KO_all_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = DP_KO_all_anno$feature, column = "GENENAME", keytype = "GENEID")
pre1_KO_all_anno$feature[is.na(pre1_KO_all_anno$feature)] <- "."
pre1_KO_all_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = pre1_KO_all_anno$feature, column = "GENENAME", keytype = "GENEID")
pre2_KO_all_anno$feature[is.na(pre2_KO_all_anno$feature)] <- "."
pre2_KO_all_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = pre2_KO_all_anno$feature, column = "GENENAME", keytype = "GENEID")
mature_KO_all_anno$feature[is.na(mature_KO_all_anno$feature)] <- "."
mature_KO_all_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = mature_KO_all_anno$feature, column = "GENENAME", keytype = "GENEID")
DP_WT_anno$feature[is.na(DP_WT_anno$feature)] <- "."
DP_WT_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = DP_WT_anno$feature, column = "GENENAME", keytype = "GENEID")
pre1_WT_anno$feature[is.na(pre1_WT_anno$feature)] <- "."
pre1_WT_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = pre1_WT_anno$feature, column = "GENENAME", keytype = "GENEID")
pre2_WT_anno$feature[is.na(pre2_WT_anno$feature)] <- "."
pre2_WT_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = pre2_WT_anno$feature, column = "GENENAME", keytype = "GENEID")
mature_WT_anno$feature[is.na(mature_WT_anno$feature)] <- "."
mature_WT_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = mature_WT_anno$feature, column = "GENENAME", keytype = "GENEID")
DP_KO_anno$feature[is.na(DP_KO_anno$feature)] <- "."
DP_KO_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = DP_KO_anno$feature, column = "GENENAME", keytype = "GENEID")
pre1_KO_anno$feature[is.na(pre1_KO_anno$feature)] <- "."
pre1_KO_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = pre1_KO_anno$feature, column = "GENENAME", keytype = "GENEID")
pre2_KO_anno$feature[is.na(pre2_KO_anno$feature)] <- "."
pre2_KO_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = pre2_KO_anno$feature, column = "GENENAME", keytype = "GENEID")
mature_KO_anno$feature[is.na(mature_KO_anno$feature)] <- "."
mature_KO_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = mature_KO_anno$feature, column = "GENENAME", keytype = "GENEID")
DP_sepcific_anno$feature[is.na(DP_sepcific_anno$feature)] <- "."
DP_sepcific_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = DP_sepcific_anno$feature, column = "GENENAME", keytype = "GENEID")
pre1a_sepcific_anno$feature[is.na(pre1a_sepcific_anno$feature)] <- "."
pre1a_sepcific_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = pre1a_sepcific_anno$feature, column = "GENENAME", keytype = "GENEID")
pre1b_sepcific_anno$feature[is.na(pre1b_sepcific_anno$feature)] <- "."
pre1b_sepcific_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = pre1b_sepcific_anno$feature, column = "GENENAME", keytype = "GENEID")
pre2a_sepcific_anno$feature[is.na(pre2a_sepcific_anno$feature)] <- "."
pre2a_sepcific_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = pre2a_sepcific_anno$feature, column = "GENENAME", keytype = "GENEID")
pre2b_sepcific_anno$feature[is.na(pre2b_sepcific_anno$feature)] <- "."
pre2b_sepcific_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = pre2b_sepcific_anno$feature, column = "GENENAME", keytype = "GENEID")
mature_sepcific_anno$feature[is.na(mature_sepcific_anno$feature)] <- "."
mature_sepcific_anno$geneName <- mapIds(EnsDb.Mmusculus.v79, keys = mature_sepcific_anno$feature, column = "GENENAME", keytype = "GENEID")

データフレームに変換してファイルの書き出し

df_DP_WT_all_anno <- as.data.frame(DP_WT_all_anno)
write.table(df_DP_WT_all_anno, "gene_name_annotation/DP_WT_all_anno.txt", sep = "\t", quote = F)
df_pre1_WT_all_anno <- as.data.frame(pre1_WT_all_anno)
write.table(df_pre1_WT_all_anno, "gene_name_annotation/pre1_WT_all_anno.txt", sep = "\t", quote = F)
df_pre2_WT_all_anno <- as.data.frame(pre2_WT_all_anno)
write.table(df_pre2_WT_all_anno, "gene_name_annotation/pre2_WT_all_anno.txt", sep = "\t", quote = F)
df_mature_WT_all_anno <- as.data.frame(mature_WT_all_anno)
write.table(df_mature_WT_all_anno, "gene_name_annotation/mature_WT_all_anno.txt", sep = "\t", quote = F)
df_DP_KO_all_anno <- as.data.frame(DP_KO_all_anno)
write.table(df_DP_KO_all_anno, "gene_name_annotation/DP_KO_all_anno.txt", sep = "\t", quote = F)
df_pre1_KO_all_anno <- as.data.frame(pre1_KO_all_anno)
write.table(df_pre1_KO_all_anno, "gene_name_annotation/pre1_KO_all_anno.txt", sep = "\t", quote = F)
df_pre2_KO_all_anno <- as.data.frame(pre2_KO_all_anno)
write.table(df_pre2_KO_all_anno, "gene_name_annotation/pre2_KO_all_anno.txt", sep = "\t", quote = F)
df_mature_KO_all_anno <- as.data.frame(mature_KO_all_anno)
write.table(df_mature_KO_all_anno, "gene_name_annotation/mature_KO_all_anno.txt", sep = "\t", quote = F)
df_DP_WT_anno <- as.data.frame(DP_WT_anno)
write.table(df_DP_WT_anno, "gene_name_annotation/DP_WT_anno.txt", sep = "\t", quote = F)
df_pre1_WT_anno <- as.data.frame(pre1_WT_anno)
write.table(df_pre1_WT_anno, "gene_name_annotation/pre1_WT_anno.txt", sep = "\t", quote = F)
df_pre2_WT_anno <- as.data.frame(pre2_WT_anno)
write.table(df_pre2_WT_anno, "gene_name_annotation/pre2_WT_anno.txt", sep = "\t", quote = F)
df_mature_WT_anno <- as.data.frame(mature_WT_anno)
write.table(df_mature_WT_anno, "gene_name_annotation/mature_WT_anno.txt", sep = "\t", quote = F)
df_DP_KO_anno <- as.data.frame(DP_KO_anno)
write.table(df_DP_KO_anno, "gene_name_annotation/DP_KO_anno.txt", sep = "\t", quote = F)
df_pre1_KO_anno <- as.data.frame(pre1_KO_anno)
write.table(df_pre1_KO_anno, "gene_name_annotation/pre1_KO_anno.txt", sep = "\t", quote = F)
df_pre2_KO_anno <- as.data.frame(pre2_KO_anno)
write.table(df_pre2_KO_anno, "gene_name_annotation/pre2_KO_anno.txt", sep = "\t", quote = F)
df_mature_KO_anno <- as.data.frame(mature_KO_anno)
write.table(df_mature_KO_anno, "gene_name_annotation/mature_KO_anno.txt", sep = "\t", quote = F)
df_DP_sepcific_anno <- as.data.frame(DP_sepcific_anno)
write.table(df_DP_sepcific_anno, "gene_name_annotation/DP_specific_anno.txt", sep = "\t", quote = F)
df_pre1a_sepcific_anno <- as.data.frame(pre1a_sepcific_anno)
write.table(df_pre1a_sepcific_anno, "gene_name_annotation/pre1a_specific_anno.txt", sep = "\t", quote = F)
df_pre1b_sepcific_anno <- as.data.frame(pre1b_sepcific_anno)
write.table(df_pre1b_sepcific_anno, "gene_name_annotation/pre1b_specific_anno.txt", sep = "\t", quote = F)
df_pre2a_sepcific_anno <- as.data.frame(pre2a_sepcific_anno)
write.table(df_pre2a_sepcific_anno, "gene_name_annotation/pre2a_specific_anno.txt", sep = "\t", quote = F)
df_pre2b_sepcific_anno <- as.data.frame(pre2b_sepcific_anno)
write.table(df_pre2b_sepcific_anno, "gene_name_annotation/pre2b_specific_anno.txt", sep = "\t", quote = F)
df_mature_sepcific_anno <- as.data.frame(mature_sepcific_anno)
write.table(df_mature_sepcific_anno, "gene_name_annotation/mature_specific_anno.txt", sep = "\t", quote = F)

特異的な遺伝子を探す
developmental stageサンプル

cut -f 20 macs/gene_name_annotation/DP_specific_anno.txt | sort -u > DP_specific_symbol_only.txt
cut -f 20 macs/gene_name_annotation/pre1a_specific_anno.txt | sort -u > pre1a_specific_symbol_only.txt
cut -f 20 macs/gene_name_annotation/pre1b_specific_anno.txt | sort -u > pre1b_specific_symbol_only.txt
cut -f 20 macs/gene_name_annotation/pre2a_specific_anno.txt | sort -u > pre2a_specific_symbol_only.txt
cut -f 20 macs/gene_name_annotation/pre2b_specific_anno.txt | sort -u > pre2b_specific_symbol_only.txt
cut -f 20 macs/gene_name_annotation/mature_specific_anno.txt | sort -u > mature_specific_symbol_only.txt

WTとKOの比較

for input in DP_KO_anno DP_WT_all_anno DP_WT_anno mature_KO_all_anno mature_KO_anno mature_WT_all_anno mature_WT_anno pre1_KO_all_anno pre1_KO_anno pre1_WT_all_anno pre1_WT_anno pre2_KO_all_anno pre2_KO_anno pre2_WT_all_anno pre2_WT_anno; do cut -f 20 macs/gene_name_annotation/${input}.txt | sort -u > ${input}_symbol_only.txt ; done

### find motifs

for input in {1..16} ; do findMotifsGenome.pl macs2/sample${input}_summits.bed mm10 ./motif/sample${input} -size 200 -mask -p 8 ; done

for input in DP_KO.bed mature_specific.bed pre2_KO.bed DP_WT.bed pre1_KO.bed pre2_WT.bed DP_specific.bed pre1_WT.bed pre2a_specific.bed mature_KO.bed pre1a_specific.bed pre2b_specific.bed mature_WT.bed pre1b_specific.bed ; do findMotifsGenome.pl macs2/specific/${input} mm10 ./motif/specific/sample${input} -size 200 -mask -p 8 ; done

https://tobiasrausch.com/courses/atac/index.html#prepare-the-input-count-matrix
https://rockefelleruniversity.github.io/RU_ATAC_Workshop.html
https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1009865#sec010