chip-seq 大麦 V2版本参考基因组

1.参考基因组

1.1 fasta文件

在/u2/new2020barley/barley_annotion 这个位置

image.png

1.2 gff文件

在/u2/new2020barley/barley_annotion/gene_annotation_morex_v2/mascher_IPK-GATERSLEBEN.DE/barley_annotation_morexv2/gene_annotation 这个位置

image.png

2.bowtie2-build 对参考基因组建立索引

bowtie2-build Hordeum_vulgare.MorexV3_pseudomolecules_assembly.dna.toplevel.fa MorexV3
##bowtie2-build 中间没有空格

image.png

3.bowtie2 进行比对

SE50利用bowtie2进行单端比对

for i in find /u2/chip_seq.test/clean/clean_data -name *V1*`;do bowtie2 -p 20 -x /u2/chip_seq.test/barley_morex.v2.annotation/Morex_V2 ${i} -S ${i}.sam;done


##sam文件变成bam文件
for i in `find /u2/chip_seq.test/morex_v2_result/2.bowtie2_result -name *.sam`;do /home/shenkuocheng/conda3/bin/samtools view -bS ${i} > ${i}.bam;done
##samtools sort
for i in `find /u2/chip_seq.test/morex_v2_result/2.bowtie2_result -name *.bam`;do /home/shenkuocheng/conda3/bin/samtools sort -@ 8 ${i} -o  ${i}.sorted.bam;done

##对bam文件index
for i in `find /u2/chip_seq.test/morex_v2_result/2.bowtie2_result -name *sorted.bam`;do samtools index ${i};done

image.png

4.去除PCR重复

sambamba markdup -t 2 -r V1-B-anti_2_1.sorted.bam V1-B-anti_2_1.markdup.bam >V1-B-anti_2_1.log 2>&1
sambamba markdup -t 2 -r V1-B-input_1.sorted.bam V1-B-input_1.markdup.bam >V1-B-input_1.log 2>&1
sambamba markdup -t 2 -r V1-B-neg_2_1.sorted.bam V1-B-neg_2_1.markdup.bam >V1-B-neg_2_1.log 2>&1
sambamba markdup -t 2 -r V1-M-anti_2_1.sorted.bam V1-M-anti_2_1.markdup.bam >V1-M-anti_2_1.log 2>&1
sambamba markdup -t 2 -r V1-M-input_1.sorted.bam V1-M-input_1.markdup.bam >V1-M-input_1.log 2>&1
sambamba markdup -t 2 -r V1-M-neg_1.sorted.bam V1-M-neg_1.markdup.bam >V1-M-neg_1.log 2>&1



V1-B-neg_2_1

5.计算基因组的有效大小

大麦基因组大小为4.2e9

#coding=utf-8
import sys
aList=[]
fa_file = sys.argv[1]
with open(fa_file,'r') as f:
    for line in f:
        line = line.strip()
        line = line.upper()
        if not line.startswith(">"):
            baseA = line.count("A")
            baseT = line.count("T")
            baseC = line.count("C")
            baseG = line.count("G")
            aList.extend([baseA, baseT, baseC, baseG])
            # print(aList)
    print("有效基因组大小：", sum(aList))

(16条消息) 计算基因组大小awk_bioinfo的博客-CSDN博客计算基因组大小

image.png

6.macs2找peak峰

macs2 callpeak -t V1-B-anti_2_1.markdup.bam -c V1-B-neg_2_1.markdup.bam -f BAM -g 4.2e9 -n V1-B-anti_neg -q 0.05 -B >V1-B-anti_neg.log 2>&1
 macs2 callpeak -t V1-B-anti_2_1.markdup.bam -c V1-B-input_1.markdup.bam -f BAM -g 4.2e9 -n V1-B-anti_input -q 0.05 -B >V1-B-anti_input.log 2>&1
macs2 callpeak -t V1-M-anti_2_1.markdup.bam -c V1-M-neg_1.markdup.bam -f BAM -g 4.2e9 -n V1-M-anti_neg -q 0.05 -B >V1-M-anti_neg.log 2>&1
macs2 callpeak -t V1-M-anti_2_1.markdup.bam -c V1-M-input_1.markdup.bam -f BAM -g 4.2e9 -n V1-M-anti_input -q 0.05 -B >V1-M-anti_input.log 2>&1

# -f 指定输入文件的格式，如SAM、BAM、BED等
# -c 对照组，可以接多个数据，空格分隔。
# -t 实验组，ChIP-seq数据。可以接多个数据，空格分隔。
# -n 输出文件名的前缀
# -g 有效基因组大小。软件有几个默认值，如hs指human
# --outdir 输出结果的存储路径
# --bdg 输出 bedgraph 格式的文件
# -q FDR阈值, 默认为0.05

MACS2 call peaks的统计学原理 - 简书 (jianshu.com)

7.HOMER查找motif

##利用conda 安装homer
# 在bioconda中搜索homer，可以看到有很多种版本
conda search homer -c bioconda

# 默认安装 homer 最新版
conda install -c bioconda homer -y # -y 表示直接安装，不询问确认与否

显示安装成功

image.png

但是后面查看有没有大麦的基因组文件的时候，需要用的configureHomer.pl这个程序，发现安装在anaconda3/share/homer 目录下的文件，修改环境配置

vim /etc/profile
export PATH="/root/anaconda3/share/homer/bin:$PATH"
##添加后source一下
source /etc/profile
##这个程序让我误删了，需要在重新下载，用wget 下载configureHomer.pl
wget http://homer.ucsd.edu/homer/configureHomer.pl

这个需要单独放在一个文件里面
/root/anaconda3/share/homer/configureHomer.pl

我用的大麦的注释是V2版本进行call peak，Ensemble中的版本是V3版本，命名规则是不一样的。这样的话这个软件用起来就比较麻烦了。

image.png

8.IGV软件展示bam文件

这个步骤应该在前面进行，自己安装的IGV一直有毛病，卸载后重新安装后又继续进行这一部分的学习
IGV软件是查看测序深度进行可视乎的软件。对chipseq数据用IGV展示的时候挑选峰的准则是 其他组的峰在对照组中没有，则为peak。
Morex_V2_scaf - Genome - Assembly - NCBI (nih.gov)

玩转基因组浏览器之tdf文件 - 云+社区 - 腾讯云 (tencent.com)
科研进展-上海天昊生物科技有限公司|上海天昊遗传分析中心 (geneskybiotech.com)