2019-04-24-ATAC-seq

for i in *_1.fastq ; do i=${i%_1.fastq*}; perl '/media/shen/disk1/biosoft/TrimGalore-0.6.0/trim_galore' --paired --gzip ${i}_1.fastq ${i}_2.fastq -o ./ ; done

bowtie2 -p 40 -t -q -N 1 -L 25 -X 2000 --no-mixed --no-discordant -x '/media/shen/disk1/jjc/reference/hg19/hg19' -1 SRR5837329_1_val_1.fq.gz -2 SRR5837329_2_val_2.fq.gz -S SRR5837329.sam ; bowtie2 -p 40 -t -q -N 1 -L 25 -X 2000 --no-mixed --no-discordant -x '/media/shen/disk1/jjc/reference/hg19/hg19' -1 SRR5837330_1_val_1.fq.gz -2 SRR5837330_2_val_2.fq.gz -S SRR5837330.sam

cd /media/shen/disk1/jjc/2019-5-7-atac-seq-mouse-mm9/no-paramater/ ; for i in *.sam; do i=${i%.sam*}; samtools sort -@ 40 -o ${i}.bam ${i}.sam ; done ; for i in *.bam; do i=${i%.bam*}; java -jar ~/下载/picard.jar MarkDuplicates REMOVE_DUPLICATES=true I=${i}.bam O=${i}.dedup.bam M=${i}.metrix.txt ; done

###去dup##

for i in *.bam; do i=${i%.bam*}; java -jar ~/下载/picard.jar MarkDuplicates REMOVE_DUPLICATES=true I=${i}.bam O=${i}.dedup.bam M=${i}.metrix.txt ; done

###去除线粒体##

samtools view -h ${i}.dedup.bam | grep -v 'chrM' | samtools view -bS -o ${i}.final.bam

### deeptools 处理###

bamCoverage -b reads.bam -o coverage.bw -bs 1000 -p 30 --normalizeUsing RPKM -e 250 

multiBigwigSummarybins-bfile1.bwfile2.bw-oresults.npz

plotCorrelation -in female_multiBigwigSummary.npz -o female_heatmap.png --corMethod spearman -p heatmap

plotCorrelation -in female_multiBigwigSummary.npz -o female_scatter.png --corMethod spearman -p scatterplot

###选择profile的 范围##

bedtools intersect -a /media/shen/disk1/jjc/reference/mm9/mm9.LINE.bed -b female_wt.macs2_peaks.narrowPeak -wa > female_wt.LINE.bed

cat female_wt.LINE.bed female_ko.LINE.bed > allfemale_line.bed

sort -k1,1 -k2,2n allfemale_line.bed > sortedallfemale_line.bed

uniq sortedallfemale_line.bed > allfemale_line.bed

computeMatrix reference-point --referencePoint TSS -S<biwigfile(s)> -R<bedfile> -p max --skipZeros -a 3000 -b 3000 -o matrix.gz 

plotProfile -m matrix.gz  -out ExampleProfile1.png (--perGroup --yMax 0.25 --yMin 0.00) --plotWidth 9 --colors red blue

###Genrich callpeak##

cd /media/shen/disk3/biosoft/Genrich-master/

./Genrich -t 1.bam,2.bam -o x.narrowPeak -j -v -e chrM

##macs2 callpeak###

macs2 callpeak --nomodel --nolambda -t wy-0331-5.dedup.bam wy-0331-6.dedup.bam wy-0331-7.dedup.bam -g mm -n female_ko --outdir female_ko/ ; macs2 callpeak --nomodel --nolambda -t wy-0331-8.dedup.bam wy-0331-9.dedup.bam wy-0331-10.dedup.bam -g mm -n male_ko --outdir male_ko/ ;

macs2 callpeak -g mm -n female_wt.macs2 -B -q 0.05 -t wy-0331-1.dedup.bam wy-0331-2.dedup.bam --nomodel --shift=25 --SPMR > callpeak.log 2>&1

bedtools intersect -a female_ko_peaks.narrowPeak -b female_ctrl_peaks.narrowPeak -v > female_koonly.narrowPeak

#macs14##

macs14 -t female_ko.bam -g mm -n female_ko.macs14 -B --nolambda --nomodel

##homer 注释##

source /etc/profile

annotatePeaks.pl male_ko.bed mm10 -cpu 20 > male_ko_anno.txt -annStats male_ko_annostat.txt

##ChIpseeker注释##

txdb <- makeTxDbFromGFF("Mus_musculus.GRCm38.94.chr.gtf",format="gtf")

male_wt <- readPeakFile("")

male_ko <- readPeakFile("")

peaks <- list(male_wt=male_wt,male_ko=male_ko)

tagMatrixList <- lapply(peaks, getTagMatrix, windows=promoter)

tagHeatmap(tagMatrixList, xlim=c(-3000, 3000), color=NULL)

peakAnnoList <- lapply(peaks, annotatePeak, TxDb=txdb,tssRegion=c(-3000, 3000), verbose=FALSE)

plotAvgProf(tagMatrixList, xlim=c(-3000, 3000))

##find overlap##

library(ChIPpeakAnno)

ol <- findOverlapsOfPeaks(male_ko,male_wt)

png('male_overlapVenn.png')

makeVennDiagram(ol)

dev.off()

###R package ATACseqQC ##

setwd("dedup/")

bamfile <- grep("dedup.bam",list.files(directory),value=TRUE)

directory <- '/media/shen/disk2/jjc/2019-4-24-atac-seq-mouse/align-results/dedup/'

bamfile.sample.ID <- gsub(".dedup.bam", "", basename(bamfile))

pdf(paste0(bamfile.sample.ID, ".fragment.size.distribution.pdf"), width =10, height=8)

fragSize <- fragSizeDist(bamFiles=bamfile, bamFiles.labels=bamfile.sample.ID)

dev.off()