这里是佳奥！让我们继续吧！

这里使用的示例文件是：antiH3K27ry-1和antiH3K27ry-2，直接Google就可以找到。

##antiH3K27ry-1
https://trace.ncbi.nlm.nih.gov/Traces/index.html?view=run_browser&acc=SRR3102823&display=metadata

##antiH3K27ry-2
https://trace.ncbi.nlm.nih.gov/Traces/index.html?view=run_browser&acc=SRR3102824&display=metadata

不过呢，这里点击FASTA/FASTQ download可以直接下fastq文件，本着从简想法，就直接在浏览器下载了。

1 质量控制与过滤

##新建qc目录
$ ls ../raw_fq/*gz| xargs fastqc -t 10 -o  ./

查看.html结果：length只有51，所以过滤选择35的话会过滤太多reads，这里选择25。

QQ截图20220823114839.png QQ截图20220823114857.png

这部分需要进行质量控制。

trim_galore软件需要依赖cutadapt。

##确保安装cutadapt
conda install -c bioconda cutadapt

##添加软件到环境变量
export PATH="$PATH:/home/kaoku/biosoft/trimgalory/TrimGalore-0.6.7"

##trim_galore软件进行过滤，这里的是单端测序
analysis_dir=/home/kaoku/project/fly/chip-seq

ls /home/kaoku/project/fly/chip-seq/raw_fq/*.gz | while read fq1;
do 
trim_galore -q 25 --phred33 --length 25 -e 0.1 --stringency 4 -o $analysis_dir/clean_fq  $fq1 & 
done

原本qc目录清空，再做一次fastqc

##进入qc目录
$ ls ../clean_fq/*gz| xargs fastqc -t 10 -o  ./

QQ截图20220823121153.png QQ截图20220823124613.png

结果来看好了不少。

对比一下质控前后变化：还是少了很多的。

##过滤前
(rnaseq) root 12:44:07 /home/kaoku/project/fly/chip-seq/raw_fq
$ ls -lh
总用量 608M
-rw-r--r-- 1 kaoku kaoku 548M  8月 22 17:26 antiH3K27ry-1.fastq.gz
-rw-r--r-- 1 kaoku kaoku  60M  8月 22 17:26 antiH3K27ry-2.fastq.gz

##过滤后
(rnaseq) root 12:43:22 /home/kaoku/project/fly/chip-seq/clean_fq
$ ls -lh
总用量 300M
-rw-r--r-- 1 root root 2.9K  8月 23 12:03 antiH3K27ry-1.fastq.gz_trimming_report.txt
-rw-r--r-- 1 root root 281M  8月 23 12:03 antiH3K27ry-1_trimmed.fq.gz
-rw-r--r-- 1 root root 2.7K  8月 23 12:01 antiH3K27ry-2.fastq.gz_trimming_report.txt
-rw-r--r-- 1 root root  19M  8月 23 12:01 antiH3K27ry-2_trimmed.fq.gz

样品多的话，就用multiqc整合数据：

##进入qc目录
multiqc ./

$ multiqc ./

  /// MultiQC 🔍 | v1.13.dev0

|           multiqc | Search path : /home/kaoku/project/fly/chip-seq/qc
|         searching | ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 4/4
|            fastqc | Found 2 reports
|           multiqc | Compressing plot data
|           multiqc | Report      : multiqc_report.html
|           multiqc | Data        : multiqc_data
|           multiqc | MultiQC complete

QQ截图20220823125210.png

注意一些软件会对python版本起冲突，一般专门建立一个python2环境以运行在默认python3环境冲突的软件（比如rnaseq rnaseq2环境之类）。最好在base环境只安装自己有把握不冲突的软件。

2 比对，去除PCR重复

2.1 对质控后的序列进行比对

##进入align目录
bowtie2_index=/home/kaoku/project/fly/refer/bowtie2-index/BDGP.

ls ../clean_fq/*gz | while read id;
do 
file=$(basename $id )
sample=${file%%.*}
echo $file $sample

bowtie2  -p 5  -x  $bowtie2_index -U  $id | samtools sort  -O bam  -@ 5 -o - > ${sample}.bam 
done 

2.2 合并.bam

##新建目录
mkdir mergeBam

cd ../align
ls *.bam| sed 's/_[0-9]_trimmed.bam//g' |sort -u | while read id; do samtools merge ../mergeBam/$id.merge.bam $id*.bam ; done

2.3 去除PCR重复

ls  *merge.bam  | while read id ; do ( samtools markdup -r $id $(basename $id ".bam").rmdup.bam & );done

##想要调试代码，加一个echo
ls  *merge.bam  | while read id ; do (echo samtools markdup -r $id $(basename $id ".bam").rmdup.bam & );done

2.4 构建.bam文件的index

ls  *.merge.bam  | xargs -i samtools index {} 

ls  *.rmdup.bam  | xargs -i samtools index {} 

3 寻找peaks

3.1 macs2 callpeak

##合并的bam文件运行macs2（去除PCR重复的样本），dm是果蝇的意思
ls  *merge.rmdup.bam |cut -d"." -f 1 |while read id;
do 
    if [ ! -s ${id}_rmdup_summits.bed ];
    then 
echo $id 
macs2 callpeak -c  Control.merge.rmdup.bam -t $id.merge.rmdup.bam -f BAM -B -g dm -n $id --outdir ../peaks  2> $id.log &  
    fi 
done  

补充：

除了conda安装，还可以使用pip install macs2，但是我没用过，好处是占用空间小，但是对于python等软件版本要求高，比如会遇到python version must >= 3.5之类的要求。

至此ChIP-Seq上游分析结束。

下一篇我们继续RNA-Seq的分析。

我们下一篇再见！