Chipseq单端测序数据分析

工作目录

mkdir -p ~/maos/chipseq
cd ~/maos/chipseq
mkdir {sra,fastq,fastqc,trim,bam,bw,bed,peak,pic}

download fastq file

#method 1: sratoolkits
cat srr.txt | while read line
do
    prefetch -O $wkd/sra $line
done
#method 2: axel
#find the links from ebi (srr)
cat link.txt | while read line
do
    axel -n 30 ftp://${line}
done

qc & trim adapter

ls fastq/*.fastq.gz |xargs fastqc -t $nThread -o fastqc
ls fastq/*fastq.gz | while read id
do
        sample=$(basename $id ".fastq.gz")
        #echo $sample
        if [ ! -e trim/${sample}_trimmed.fq.gz ];then
                echo "trim_galore --gzip -o trim $id"
        fi
done

bowtie2 mapping

nThread=32
bt2_hg38=~/reference.and.annotations/Homo_sapiens/UCSC/hg38/Sequence/Bowtie2Index/hg38
ls trim/*.fq.gz | while read file
do
    id=$(basename $file _trimmed.fq.gz)
    echo $id
    #sample=${id%%.*}
    #echo $sample
    bowtie2 -p $nThread -k 1 -x $bt2_hg38 -U $file | samtools sort -O bam -o bam/${id}.bam
done

samtools merge

samtools merge -@ $nThread WT.bam 1.bam 2.bam
samtools merge -@ $nThread KO.bam 3.bam 4.bam

bam2bw

# no input
ls bam/*bam |while read file
do
    samtools index $file $file.bai
    id=$(basename $file .bam)
    echo trans $file to ${id}.bw
    bamCoverage -b $file -o bw/${id}.bw -p $nThread
done
# with input
ls bam/*bam |while read file
do
    samtools index -@ $nThread $file $file.bai
done
bamCompare -p $nThread --bamfile1 bam/SMARCA4_H.bam --bamfile2 bam/input_H.bam --outFileName bw/SMARCA4_H.bw
bamCompare -p $nThread --bamfile1 bam/SMARCA4_P.bam --bamfile2 bam/input_P.bam --outFileName bw/SMARCA4_P.bw
bamCompare -p $nThread --bamfile1 bam/SMARCC1_H.bam --bamfile2 bam/input_H.bam --outFileName bw/SMARCC1_H.bw
bamCompare -p $nThread --bamfile1 bam/SMARCC1_P.bam --bamfile2 bam/input_P.bam --outFileName bw/SMARCC1_P.bw

macs2 callpeak

ls bam/*.bam | while read file
do
    id=$(basename $file .bam)
    macs2 callpeak -t $file -m 10 30 -p 1e-5 -f BAMPE -g hs -n ${id} --outdir peak 2>${id}.masc2.log
done