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)
网友评论