热点新闻
ATAC-bulk rna的染色质可及性分析流程
2024-11-09 21:59  浏览:504  搜索引擎搜索“微商筹货网”
温馨提示:信息一旦丢失不一定找得到,请务必收藏信息以备急用!本站所有信息均是注册会员发布如遇到侵权请联系文章中的联系方式或客服删除!
联系我时,请说明是在微商筹货网看到的信息,谢谢。
展会发布 展会网站大全 报名观展合作 软文发布

00.背景

基于常规转录组(bulk rna)的数据进行染色质可及性(ATAC)的常规分析流程。代码仅供参考。

01.质控

和别的组学一样,拿到原始fastq数据的第一步时进行质控处理,这里可以选择的软件是trim_galore或者trimmomatic。首先需要去相关软件的官网或者是其他教程介绍中下载相应的软件。对于trim_galore,参考代码如下:

nohup trim_galore -q 25 --phred33 --length $minlen -e 0.1 --stringency 4 --paired -o ./02.clean $fq1 $fq2 & #'$minlen质控时候的最小长度 #'$fq1 左端序列 #'$fq2 右端序列

02.比对

将序列与对应的参考基因组进行比对,为序列注释信息。这里可以选择软件bowtie2和samtools,参考代码如下:

if [ ! -d "03.align" ];then mkdir 03.align samtools=软件路径 cut -f1 sample.list|rush -j 8 \ "bowtie2软件路径 \ -1 ./02.clean/{}_1.clean.fq.gz -2 ./02.clean/{}_2.clean.fq.gz \ -x $bt2_genome|$samtools sort -@ 5 -O bam -o ./03.align/{}.raw.bam " #比对完马上生成索引 for i in `cut -f1 sample.list`;do $samtools index -@ 5 -b ./03.align/$i.raw.bam $samtools flagstat ./03.align/$i.raw.bam >./03.align/$i.raw.stat done else echo "03.align already exists." fi #'sample.list是一个信息列表,主要为了方便传参并行分析 #'相应的参数可以去官网或者相关教程查看

03.去重

ATAC数据的特点是需要去除PCR重复,可以使用的软件是picard。同时,线粒体和低质量的序列也需要过滤。
去重后统计每个样本的序列信息:

##mapping quality>30;注意线粒体 for i in `cut -f1 sample.list`;do samtools view -h -f 2 -q 30 ./03.align/$i.rmdup.bam \ | grep -v -e "chrM" \ | samtools sort -O bam -@ 10 -o - > ./03.align/$i.last.bam samtools index ./03.align/$i.last.bam samtools flagstat ./03.align/$i.last.bam > ./03.align/$i.last.stat done

使用python脚本处理质控信息表格:

import sys import os import pandas as pd sample_list = sys.argv[1] def get_FQinfo(sample_list): with open(f"{sample_list}", 'r') as f: f = f.read().split('\n') while '' in f: f.remove('') sample_id=[] final_reads=[] duplicate=[] read1_mapped=[] read2_mapped=[] for line in f: line = line.split('\t') sample_name = line[0] sample_id.append(sample_name) with open(f"./03.align/{sample_name}.last.stat", 'r') as f1: f1 = f1.read() f1 = f1.split('\n') final_reads.append(str.split(f1[0]," +")[0]) read1_mapped.append(str.split(f1[6]," +")[0]) read2_mapped.append(str.split(f1[7]," +")[0]) duplicate.append(str.split(f1[3]," +")[0]) qc_table=pd.Dataframe() qc_table["Sample ID"]=sample_id qc_table["Read1 Mapped"]=read1_mapped qc_table["Read2 Mapped"]=read2_mapped qc_table["Final Reads"]=final_reads qc_table["Duplicates"]=duplicate qc_table.to_csv("./03.align/align_qc_stat.csv") get_FQinfo(sample_list)

04.peaks分析

主要目的是找到序列中的peaks,用到的软件是macs2。参考代码如下:

cut -f1 sample.list|rush -j 8 \ "macs2 callpeak -t ./04.peaks/{}.last.bed -g mm --nomodel --shift -100 --extsize 200 -n {} --outdir ./04.peaks/"

05.TSS/TES可视化

需要使用deeptools将比对的bam文件转换为bw文件。
然后使用computeMatrix进行分析并使用plotHeatmap和plotProfile进行可视化。
参考代码如下:

cut -f1 sample.list|rush -j 8 \ " mkdir -p ./04.peaks/TSS/{} computeMatrix reference-point --referencePoint TSS -p 15 \ -b 2000 -a 2000 \ -R refseq.bed \ -S ./04.peaks/{}.last.bw \ -o ./04.peaks/TSS/{}/TSS.gz \ --skipZeros \ --outFileSortedRegions ./04.peaks/TSS/{}/Heatmap1sortedRegions.bed "

for i in `cut -f1 sample.list`;do plotHeatmap -m ./04.peaks/TSS/$i/TSS.gz -out ./04.peaks/TSS/$i/$i_TSS.Heatmap.pdf --plotFileFormat pdf plotProfile -m ./04.peaks/TSS/$i/TSS.gz -out ./04.peaks/TSS/$i/$i_TSS.Profile.pdf --plotFileFormat pdf --perGroup done

06.样本相关性

主要使用deeptools进行样本相关性分析和可视化。

07.peaks结构注释

ChIPseeker注释peaks,仅仅是结构注释:启动子等。这里需要用到R脚本。

07.peaks motif注释

注释motif,使用homer软件。注意物种。

perl ~/miniconda3/envs/atac/share/homer-4.9.1-6/configureHomer.pl -install mm10

08.peaks差异分析

差异peaks,使用DiffBind软件。

09.footprint分析

使用hintatac软件,注意物种。

发布人:d01b****    IP:124.223.189***     举报/删稿
展会推荐
让朕来说2句
评论
收藏
点赞
转发