从SRA数据下载开始学习ATACseq数据分析

2023-10-20 13:20
文章标签 数据 学习 下载 数据分析 sra atacseq

本文主要是介绍从SRA数据下载开始学习ATACseq数据分析，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

1.打开NCBI GEO数据库找到你需要下载的数据页，打开SRA Run Selector，找到需要下载的原始数据的SRR编号。

 

2.创建一个txt文件并输入需要下载的SRR编号（好像SRA页面可以直接生成？待研究）

vim sra.txt

3.prefetch命令进行数据下载（很慢，所以要放在后台下载，并且关闭shell窗口不影响）

nohup prefetch  --option-file ./sra.txt &

下载好之后会在本地出现以SRR编号命名的文件夹，需要进入其中的每一个文件进行解压 

fasterq-dump  -p -e 24 --split-3 *.sra

可以看到每个.sra文件解压成了两个.fastq文件，分别是_1.fastq和 _2.fastq。

每个fastq文件内容格式有4行：第1行主要储存序列测序时的坐标等信息；第2行储存的是序列信息，正常情况都是用ATCG四个字母表示，但是当测序仪无法准确分辨该位置的序列信息时，会以N来代指此处的序列信息；第3行以“+”开始，可以存储一些信息，但是目前这一行都是空的；第4行存储的就是第2行每一个碱基的测序质量信息，其中的每一个符号所对应计算机的ASCII值是经过换算的phred值，而phred值等于33-10*logP，这里的P代表该位置测序发生错误的概率，简单来说，如果某个位置测得的序列十分可信，那么意味着该位置发生错误的概率极小，所以phred值就很大，即该值越大，说明测序的质量越好。 

将所有生成的双端.fastq文件整合到一个文件夹中

4.fastqQC+multiQC

less /share/home/wuqian/job/atac_test/atac_test.txt |while  read id;
do 
fastqc ${id}_1.fastq -o /share/home/wuqian/job/atac_test/fastqc
fastqc ${id}_2.fastq -o /share/home/wuqian/job/atac_test/fastqc
done

获得了.html的QC文件，此时可以吧这些文件拷贝到本地进行查看。

# 横轴是1 - 97 bp；纵轴是百分比
# 图中四条线代表A T C G在每个位置平均含量
# 理论上来说，A和T应该相等，G和C应该相等，但是一般测序的时候，刚开始测序仪状态不稳定，很可能出现上图的情况。像这种情况，即使测序的得分很高，也需要cut开始部分的序列信息，那根据这张图，我准备左边cut 18bp。

根据质控报告对每个fastq文件进行质控 

##conda 安装fastp
conda install -c bioconda fastp##创建一个质控文件夹
mkdir fastp##这时我们需要写一个循环来批量质控
##逻辑是，先读取我们之前的下载SRA的文件名，对每个名字的_1的fastq和_2的fastq利用fastp进行质控
## -f 从5‘剪切碱基数目; -t 从3‘剪切碱基数目
nohup
less /share/home/wuqian/job/atac_test/atac_test.txt |while  read id;
do 
fastp -i ${id}_1.fastq -o ./fastp/${id}_1.fastq -I ${id}_2.fastq -O ./fastp/${id}_2.fastq -f 18 -t 0 -L
done &##进入fastp文件夹，对处理过的数据利用fastqc重新进行质量分析
cd /share/home/wuqian/job/atac_test/fastp
mkdir fastqc_fastp##这时我们需要写一个循环来批量检查数据质量
##逻辑是，先读取我们之前的下载SRA的文件名，对每个名字的_1的fastq和_2的fastq进行检查
##-o的意思是fastqc输出的结果都输出到/home/xuyu/atac_ara/atac_data/fastp/fastqc_fastp文件夹内
nohup 
less /share/home/wuqian/job/atac_test/atac_test.txt |while  read id;
do 
fastqc ${id}_1.fastq -o /share/home/wuqian/job/atac_test/fastp/fastqc_fastp
fastqc ${id}_2.fastq -o /share/home/wuqian/job/atac_test/fastp/fastqc_fastp
done &##将.html文件重新拷贝到本地查看QC结果

 5.reads mapping 到 genome

对测序获得的reads进行mapping，才能知道基因组上哪些区域富集比较多的reads，也就是开放的染色质区域。

mkdir genome
cd genome##下载基因组注释文件
nohup wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/635/GCF_000001635.27_GRCm39/GCF_000001635.27_GRCm39_genomic.fna.gz &gzip -d GCF_000001635.27_GRCm39/GCF_000001635.27_GRCm39_genomic.fna.gz##建立基因组索引
推荐直接在bowtie2官网下载对应的index文件，而不要尝试自己下载，下载完之后一定要知道索引文件放在什么位置#将去掉接头后的fastq文件比对会基因组，同时将sam文件转bam文件（因为sam文件太大，占用太多储存空间且不利于做大样本分析，因此需要转换为仅20-30%的bam文件），在bam文件进行sorted（按照基因组位置，这样有利于在IGV中查看变异）
##These parameters ensured that fragments up to 2 kb were allowed to align (-X2000) and that only unique aligning reads were collected (-m1).#!/bash/binbowtie2 -p 5 --very-sensitive -X 2000 -x /share/home/wuqian/job/atac_test/cleandata/index/mm9 -1 test_1.fastq -2 test_2.fastq > test.sam 
提交bash脚本samtools view -bS test.sam| samtools sort -O bam -@ 5 -o - > test.raw.sorted.bam

6.对生成的bam文件进行数据清洗

##去除线粒体基因组序列，因为在线粒体基因组中没有感兴趣的ATAC-seq peaks
##只保留两条reads比对到同一条染色体（proper paired），还有高质量的比对结果（MAPQ>=30）nohup less /share/home/wuqian/job/atac_test/atac_test.txt |while  read id; 
do
samtools view -@ 5 -f 2 -q 30 -h ./sam_bam/${id}.bam | grep -v chrM | samtools sort -@ 4  -O bam -o ./sam_bam/${id}.rmChrM.bam
done >rmChrM.log 2>&1 &##去除PCR重复序列
##markdup -r意思是把重复的都删去 -t 4 意思是使用4线程处理  nohup less /share/home/wuqian/job/atac_test/atac_test.txt |while read id; 
do
sambamba markdup -r -t 4 ${id}.rmChrM.bam ${id}_rmchrM_rmdup.bam
done >rmdup.log 2>&1 &##picard获取insert metric文件和条形图（shell脚本执行）
for i in `ls *.raw.sorted.bam`;
dosampleID="${i%.*}"echo "${sampleID}"java -Xmx32g -jar /share/home/wuqian/job/atac_test/picard.jar CollectInsertSizeMetrics I=$i OUTPUT=${sampleID}_insertMetrics H=${sampleID}_histo.pdf
doneecho $?##Tn5偏移同时将清洗之后的bam文件转为bed文件（shell脚本执行）
for i in `ls *.rmchrM.rmdup.bam`;
dosampleID="${i%.*}"echo "${sampleID}"bedtools bamtobed -i  $i | awk 'BEGIN {OFS = "\t"} ; {if ($6 == "+") print $1, $2 + 5, $3 + 5, $4, $5, $6; else print $1, $2 - 4, $3 - 4, $4, $5, $6}' > ${sampleID}.bed
doneecho $?
echo 'You have bed files!'##将ENCODE黑名单中的基因去除
for i in `ls *.bed`;
dosampleID="${i%.*}"echo "${sampleID}"bedtools subtract -a $i -b mm9.blacklist.bed > ${sampleID}_debl.bed
doneecho $?
echo 'You have bed files without blacklist regions!'

7.Peak Calling

##???对比对后的bam文件进行测序文库复杂度检验
##statc_qc.sh脚本
less /share/home/wuqian/job/atac_test/atac_test.txt | while read id;  
do bedtools bamtobed -bedpe -i ${id}_rmchrM_rmdup.bam | \awk 'BEGIN{OFS="\t"}{print $1,$2,$4,$6,$9,$10}' | sort | uniq -c | \awk 'BEGIN{mt=0;m0=0;m1=0;m2=0} ($1==1){m1=m1+1} ($1==2){m2=m2+1} {m0=m0+1} {mt=mt+$1} END{m1_m2=-1.0; if(m2>0) m1_m2=m1/m2;printf "%d\t%d\t%d\t%d\t%f\t%f\t%f\n",mt,m0,m1,m2,m0/mt,m1/m0,m1_m2}' > ${id%%.*}.nodup.pbc.qc;
donenohup bash stat_qc.sh &##输出两个文件
SRR8528251.nodup.pbc.qc
SRR8528255.nodup.pbc.qccat一下查看NRF、PBC1、PBC2值，均符合
标准为：NRF>0.9, PBC1>0.9, and PBC2>10.##将bam文件转为bed文件
nohup less /share/home/wuqian/job/atac_test/atac_test.txt |while read id; do bedtools bamtobed -i ${id}_rmchrM_rmdup.bam > ${id%%.*}.raw.bed ;done >bed.log 2>&1 &
##对于chip_seq的数据分析，拿到bam文件之后直接peak calling就可以了，对于ATAC_seq而言，一定要偏移之后再进行peak calling
##在peak calling之前首先需要将Tn5转座酶插入位点进行偏移
##对于正链上的reads需要向右偏移4bp, 比对的起始位置加4，对于负链上的reads, 则向左偏移5bp, 比对的起始位置减5bp。##macs2 peak callingnohup ls *.last.bed | while read id; do (macs2 callpeak -t $id -f BED -n "${id%%.*}" -g mm --shift -100 --extsize 200 --nomodel) ;done &##macs2 callpeak会产生3个文件：NAME_peaks.xls, NAME_peaks.narrowPeak, NAME_summits.bed，其中最有用的文件是NAME_peaks.narrowPeak 是一个纯文本BED文件，列出了每个called peak的基因组坐标，以及各种统计数据(fold-change，p值和q值等)。

8.下游可视化

（1）IGV可视化比对结果，需要bam文件转为bw文

##bed文件生成bam文件，sort后再生成bw文件
bedtools bedtobam -i test.bed -g mm9.chrom.sizes > test.last.bam
nohup samtools sort  test.last.bam -@ 16  -O BAM -o test.final.bam >nohup.log 2>&1 &
samtools index test.final.bam##生成normalized(CPM)之后的bw文件
ls *.bam |while read id ; do
nohup bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.bw &
done

（2）比较不同的peak文件

 最常使用的工具是BEDTools，比如不同生物学重复之间可以发现共同的染色质区域：bedtools intersect ；实验组和对照组比较可以发现具有差异的染色质区域：bedtools subtract

（3）peak注释

CHIPseeker R包

（4）TF motif enrichment test

HOMER 

It takes a peak file as input and checks for the enrichment of both known sequence motifs and de novo motifs

##UCSC官网下载对应的基因组注释文件，一定要是从一开始处理fastq文件开始就对应的注释文件，这里错误地使用mm39代替GRCm39，导致浪费了很多时间参数：findMotifsGenome.pl <peak/bed file> <genome file> <output directory>nohup findMotifsGenome.pl SRR8528251.bed /share/home/wuqian/job/atac_test/genome/GCF_000001635.27_GRCm39_genomic.fna motif_enrichment -preparse >motif.log 2>&1 &  nohup findMotifsGenome.pl SRR8528255.bed /share/home/wuqian/job/atac_test/genome/GCF_000001635.27_GRCm39_genomic.fna motif_enrichment -preparse >motif.log 2>&1 &

参考：https://github.com/harvardinformatics/ATAC-seq.git

明码标价之ATACseq｜生信菜鸟团

「与国同庆，万字长文」ATAC-seq实战教程：从SRA数据下载到高分辨率论文主图绘制 - 徐寅生的文章 - 知乎 https://zhuanlan.zhihu.com/p/415718382 

这篇关于从SRA数据下载开始学习ATACseq数据分析的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/247347。

相关文章

C#实现将XML数据自动化地写入Excel文件

C#实现将XML数据自动化地写入Excel文件

《C#实现将XML数据自动化地写入Excel文件》在现代企业级应用中,数据处理与报表生成是核心环节,本文将深入探讨如何利用C#和一款优秀的库,将XML数据自动化地写入Excel文件,有需要的小伙伴可以... 目录理解XML数据结构与Excel的对应关系引入高效工具：使用Spire.XLS for .NETC
阅读更多...
MySQL数据目录迁移的完整过程

MySQL数据目录迁移的完整过程

《MySQL数据目录迁移的完整过程》文章详细介绍了将MySQL数据目录迁移到新硬盘的整个过程,包括新硬盘挂载、创建新的数据目录、迁移数据（推荐使用两遍rsync方案）、修改MySQL配置文件和重启验证... 目录1，新硬盘挂载（如果有的话）2，创建新的 mysql 数据目录3，迁移 MySQL 数据（推荐两
阅读更多...
Python数据验证神器Pydantic库的使用和实践中的避坑指南

Python数据验证神器Pydantic库的使用和实践中的避坑指南

《Python数据验证神器Pydantic库的使用和实践中的避坑指南》Pydantic是一个用于数据验证和设置的库,可以显著简化API接口开发,文章通过一个实际案例,展示了Pydantic如何在生产环... 目录1️⃣ 崩溃时刻：当你的API接口又双叒崩了！2️⃣ 神兵天降：3行代码解决验证难题3️⃣ 深度
阅读更多...
MySQL快速复制一张表的四种核心方法(包括表结构和数据)

MySQL快速复制一张表的四种核心方法(包括表结构和数据)

《MySQL快速复制一张表的四种核心方法(包括表结构和数据)》本文详细介绍了四种复制MySQL表（结构+数据）的方法,并对每种方法进行了对比分析,适用于不同场景和数据量的复制需求,特别是针对超大表（1... 目录一、mysql 复制表（结构+数据）的 4 种核心方法（面试结构化回答）方法 1：CREATE
阅读更多...
详解C++ 存储二进制数据容器的几种方法

详解C++ 存储二进制数据容器的几种方法

《详解C++存储二进制数据容器的几种方法》本文主要介绍了详解C++存储二进制数据容器,包括std::vector、std::array、std::string、std::bitset和std::ve... 目录1.std::vector<uint8_t>（最常用）特点：适用场景：示例：2.std::arra
阅读更多...
SpringBoot返回文件让前端下载的几种方式

SpringBoot返回文件让前端下载的几种方式

《SpringBoot返回文件让前端下载的几种方式》文章介绍了开发中文件下载的两种常见解决方案,并详细描述了通过后端进行下载的原理和步骤,包括一次性读取到内存和分块写入响应输出流两种方法,此外,还提供... 目录01 背景02 一次性读取到内存，通过响应输出流输出到前端02 将文件流通过循环写入到响应输出流
阅读更多...
SpringBoot整合Apache Spark实现一个简单的数据分析功能

SpringBoot整合Apache Spark实现一个简单的数据分析功能

《SpringBoot整合ApacheSpark实现一个简单的数据分析功能》ApacheSpark是一个开源的大数据处理框架,它提供了丰富的功能和API,用于分布式数据处理、数据分析和机器学习等任务... 目录第一步、添加android依赖第二步、编写配置类第三步、编写控制类启动项目并测试总结ApacheS
阅读更多...
MySQL中的DELETE删除数据及注意事项

MySQL中的DELETE删除数据及注意事项

《MySQL中的DELETE删除数据及注意事项》MySQL的DELETE语句是数据库操作中不可或缺的一部分,通过合理使用索引、批量删除、避免全表删除、使用TRUNCATE、使用ORDERBY和LIMI... 目录1. 基本语法单表删除2. 高级用法使用子查询删除删除多表3. 性能优化策略使用索引批量删除避免
阅读更多...
MySQL 数据库进阶之SQL 数据操作与子查询操作大全

MySQL 数据库进阶之SQL 数据操作与子查询操作大全

《MySQL数据库进阶之SQL数据操作与子查询操作大全》本文详细介绍了SQL中的子查询、数据添加（INSERT）、数据修改（UPDATE）和数据删除（DELETE、TRUNCATE、DROP）操作... 目录一、子查询：嵌套在查询中的查询1.1 子查询的基本语法1.2 子查询的实战示例二、数据添加：INSE
阅读更多...
Linux服务器数据盘移除并重新挂载的全过程

Linux服务器数据盘移除并重新挂载的全过程

《Linux服务器数据盘移除并重新挂载的全过程》：本文主要介绍在Linux服务器上移除并重新挂载数据盘的整个过程,分为三大步：卸载文件系统、分离磁盘和重新挂载,每一步都有详细的步骤和注意事项,确保... 目录引言第一步：卸载文件系统第二步：分离磁盘第三步：重新挂载引言在 linux 服务器上移除并重新挂p
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2