ATACAseq 分析流程

参考网址：https://mp.weixin.qq.com/s/XQITfdnS6Zn6iTBtuladPw
https://github.com/DawnEve/txtBlog.py txtBlog使用介绍
http://blog.applymed.cn/index.php?k=R&id=1_0#16 WJL笔记

========================================
linux环境及软件安装
----------------------------------------
conda的下载与使用：
参考：https://mp.weixin.qq.com/s?src=11&timestamp=1622351850&ver=3099&signature=IPtzc8n2H6Xmvj0gDdoFzWHGn0hud8HMas6DtjxGXeN9gE1MKiCsA81e7aIqhp-Bl6x0ygXBc2a4WO0p0UfUmDKyihYQ-15eo0k1HAo-uETQF-uySgzRYbbSYVvMdw6y&new=1
注意：在家目录下的data目录下创建一个conda目录，在conda目录下下载conda最新版本。安装路径选择conda目录下的一个不存在目录中，否则会报错。不要将conda放在家目录下，因为家目录占的内存是服务器的，而家目录下的data目录占的是自己硬盘的内存。在重复地下载安装和卸载conda时，要把前一个conda的.condarc文件删除，不然再重新下载conda时，会提示配置文件错误（can't load configuration file, yaml格式错误等）。

step1:wget-c+网址。在data目录下用wget命令下载安装conda,安装的最后一如果选择yes则默认安装在家目录下，否则自己手动输入路径：/home/huxl/data/conda/install(注意：自己在data目录下创建了conda目录，但是未创建install目录。因为安装在/home/huxl/data/conda/目录下会报错)

step2: source .bashrc 更新一下

step3: 安装的过程中如果默认添加环境变量，则不用手动添加。如果没有，则需要手动添加环境变量。
vim .bash rc
export PATH=pwd/anaconda/bin:$PATH
退出这个文件后source .bashrc更新一下即可
step4: 找到.condarc配置文件，配置镜像。一般.condarc文件和conda的安装位置同属一个目录下。如：在家目录下安装conda,则.condarc文件也位于家目录下。
vim .condarc
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes
退出这个文件后，source .condarc更新一下即可正常使用conda

step5: conda -h可以查看conda的使用方法。

========================================
一：下载作者的数据并将sra文件转换成fastq文件（在fast_dump命令进行格式转换过程中已经把接头去掉了）
----------------------------------------

step1: 在文章中找到数据的GEO号，进入网站https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66581， 输入GEO号，下载SraRunTable.txt和 SRR_Acc_List.txt到本地，然后通过Filezila上传到服务器自己存放数据的目录。

step2：mkdir {sra,raw,clean,align,peaks,motif,qc}创建这几个目录。

step3: cd 进入sra目录，用prefetch + SRR号直接下载4个测试数据。(好像不需要知道是单端还是双端，都可以用step4的命令转换。如不需要，可直接忽略以下四个方法直接进入step 4)
如果不知道文库是single end还是paired end，如何转换成fastq？参考：https://www.cnblogs.com/RyannBio/p/9582922.html
方法一：一个简单方法是在SRA Run Browser上输入SRR编号查询: (http://trace.ncbi.nlm.nih.gov/Traces/sra/)，"Browse" -> "Run Browser" -> then input your ID，LAYOUT会标明是SINGLE还是PAIRED，也可以在Reads结果页面查看，如果展示1条read就是single end，2条reads就是paired end。

方法二：可以先用sra-stat对sra文件进行统计，从统计结果里可以知道是单末端还是双末端。

方法三：直接跑fastq-dump，用参数--split-files就自可以自动分辨出来，结果是一个文件就是single end，两个文件就是paired end。

方法四：fastq-dump的--split-spot可以用来分辨
srr="SRR3184279"
numLines=$(fastq-dump -X 1 -Z --split-spot $srr | wc -l)
if [ $numLines -eq 4 ]
then
echo "$srr is single-end"
else
echo "$srr is paired-end"
fi

step4: 将所有的下载的SRR文件合并到一个目录中： mv /home/huxl/data/datasource/shengxinpractise/sra/SRR* sra/
将某个sra文件转换成fastq文件：fastq-dump --gzip --split-3 -A SRR2927015 #跑完fast-dump已经把接头去掉了。
开多个窗口同时转换。
这是处理少量文件的方法，如果处理多个文件，建议使用shell的多文件处理方法。了解更多fastq-dump用法，参考网址：https://www.cnblogs.com/fhn7/articles/12355025.html
处理多个文件，转换为fastq格式：fastq-dump --gzip --split-3 -A $i${i}.sra && echo"** ${i}.sra to fastq done **"
--split-3 : 将双端测序分为两份,放在不同的文件,但是对于一方有而一方没有的reads会单独放在一个文件夹里；
&&：连接两个命令，前面的命令运行完毕且正确执行，后面的命令才会运行；||连接两个命令，前一个命令运行完毕且为错误，则开始执行后面的命令。|为管道符命令，后面的命令必须要能接收前一个命令的输出数据才可以。如：less, more,head,tail才可以; 而ls,cp,mv等就不行。
${}:用于区分变量的边界，明确告诉程序要取哪个变量的值。个位数的，可直接使用数字，但两位数以上，则必须使用 {} 符号来括住。如：$9, ${10}。
cho"** ${i}.sra to fastq done **“:这个命令表示如果上一个命令执行完毕，将输出已经转换完成这个提示。

========================================
二：质控（使用trim_galore命令去除低质量数据）
----------------------------------------
step1: 下载trim-galore软件

下载trimgalore软件：在下载这个软件前需要先下载fastqc软件和cutadaptor软件。
FastQC下载：
下载安装包：wget -c https://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.9.zip
解压：unzip fastqc_v0.11.9.zip
将FastQC的安装路径添加bashrc文件：export PATH=/home/huxl/data/tools/FastQC:$PATH
测试安装是否成功：fastqc --help
CutAdaptor下载：参考官网下载方式：https://cutadapt.readthedocs.io/en/stable/
（有关cutadaptor的用法原理请见：https://www.cnblogs.com/xudongliang/p/6404958.html）
使用pip一键下载并安装： python3 -m pip install --user --upgrade cutadapt
测试是否已经安装完成：cutadapt --version

Trim_Galore下载：参考官网下载方式：https://github.com/FelixKrueger/TrimGalore
使用官网给的help文件中的下载方式：curl -fsSL https://github.com/FelixKrueger/TrimGalore/archive/0.6.6.tar.gz -o trim_galore.tar.gz
解压：tar xvzf trim_galore.tar.gz
添加路径：export PATH=/home/huxl/data/tools/trim_galore.tar.gz:$PATH
测试是否安装成功：trim_galore --help
查看fastq文件里面到底包含了哪些东西？
$zcat SRR2927016_1.fastq.gz|head -n 10 #zcat和zless一样都可以查看压缩文件，head命令查看前几行， -n参数指定head查看的行数，后面接10代表查看前10行。

@SRR2927015.1 1 length=101
NGCTACACCTTGACCTAACGTTTTTATGTTTGATTCTTTTGCTTACTTTAATACCTTTTTAGGGTTTGCTGAAGATGGCGGTATATAGGCTGAATTAGCAA
+SRR2927015.1 1 length=101
#1:BDBDDFHHHFIIIGIIIAFGHIGGIIDGDHHIIGGGGIIIIIIICBFBGCEHHIIIGEHHG@@@ADEHIDAEEHFHEC3;?CDD>ACCBBCCA@>AC@
@SRR2927015.2 2 length=101
AGTCTACCCACCTCTAGCCGGAAATCTAGCCCATGCAGGAGCATCAGTAGACCTAACAATTTTCTCCCTTCATTTAGCTGGAGTGTCATCTATTTTAGGTG
+SRR2927015.2 2 length=101
B@BFDBDFHGFHHGGIIIGE>HIGJIIIHIGEIIGIIJ@BHGIJIGIF9BGGIGIJIFDIJIHHFHGGFFFFFEDECEEECDDCB@>@CCDDEEDDFDD:A
@SRR2927015.3 3 length=101
NAGTAAACTGATAGGCTAGATGTTGCTAAAATAAATAAAATCCCTAGGTTTAAATTAATTAATGGGTGTGGTATTGGTAGGGGAACTCATAGACTTAATGC
结果解读：
第一行：样本名称（第几条read），该条read长度。SRR2927015.1表示这个样本的第一条Read,经过zcat SRR2927016_1.fastq.gz|tail -n 10命令查看这个样本共有134330条read.
第二行：read序列，以ATCGN表示，N表示荧光信号干扰无法判断到底是哪个碱基。
第三行：以“+”开始，可以储存一些附加信息，一般是空的。 第四行：储存的是质量信息，与第2行的碱基序列是一一对应的，其中的每一个符号对应的ASCII值成为phred值，可以简单理解为对应位置碱基的质量值，越大说明测序的质量越好。不同的版本对应的不同。
有关碱基错误率及Q值计算方法请参考：https://blog.csdn.net/weixin_39816448/article/details/111687753

step2: 测序数据的过滤

在包含了fastq.gz的文件夹里创建config.raw文件，第一列随意填，充数用。第二列为所有fastq1文件的路径，第三列为所有fastq2文件的路径。
1 1 /home/huxl/data/datasource/shengxinpractise/sra/sra/SRR2927015_1.fastq.gz /home/huxl/data/datasource/shengxinpractise/sra/sra/SRR2927015_2. fastq.gz
2 2 /home/huxl/data/datasource/shengxinpractise/sra/sra/SRR2927016_1.fastq.gz /home/huxl/data/datasource/shengxinpractise/sra/sra/SRR2927016_2. fastq.gz
3 3 /home/huxl/data/datasource/shengxinpractise/sra/sra/SRR2927018_1.fastq.gz /home/huxl/data/datasource/shengxinpractise/sra/sra/SRR2927018_2. fastq.gz
4 4 /home/huxl/data/datasource/shengxinpractise/sra/sra/SRR3545580_1.fastq.gz /home/huxl/data/datasource/shengxinpractise/sra/sra/SRR3545580_2. fastq.gz
创建一个名为trim_galore.sh的脚本文件，使用trim-galore命令对低质量数据进行过滤。脚本文件内容如下：
mkdir -p clean #创建一个clean文件夹用作输出文件夹
cat config.raw|while read id # 和do echo $id 这条命令连用，用于查看并输出 config.raw文件
do echo $id
arr=($id) #将id作为变量赋值给arr
fq2=${arr[2]} #将arr这个列表的第三个元素，也就id的fastq2文件路径赋值给变量fq2
fq1=${arr[1]} #将arr这个列表的第二个元素，也就id的fastq1文件路径赋值给变量fq1
sample=${arr[0]}#将arr这个列表的第一个元素，也就id的1，2，3，4赋值给变sample，作为样本的名字
nohup trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paied -o clean $fq1 $fq2 & #使用trim_galore命令进行数据低质量过滤处理，-o后面接输出文件名，输入文件是变量fq1和变量fq2,输出和输入文件都在当前文件夹下，所以不用给路径。
done #写完脚本后必须写结束语句done
ps -ef |grep trim_galore #ps命令将某个进程显示出来，-f表示全格式输出
运行脚本文件：
sh trim_galore.sh
运行完成后在clean文件夹下生成
========================================
三：数据质量检测（分别对trim_galore过滤前和过滤后的数据生成网页查看）
----------------------------------------

对原始数据进行fastqc数据整合并生成网页，输出到raw文件夹：
cd ~/data/datasource/shengxinpractise/sra/sra #进入存放原始数据fastq.gz格式的文件夹下.
$mkdir -p raw #在当前目录下创建raw文件夹
$fastqc -t 5 ./*gz -o raw/
#使用fastqc命令，-t 5参数设定同时最多处理5个文件；./*gz表示输入文件为当前目录下带有gz结尾的所有文件；-o raw/表示输出文件夹为当前目录下的raw文件夹。
$multiqc ./*zip #对fastqc生成的zip文件整合并生成一个整合后的网页multiqc_report.html。

对trim_galore过滤后的数据整合并生成网页，输出到clean文件夹：
cd ~/data/datasource/shengxinpractise/sra/sra/ #进入存放trim_galore过滤后数据的文件夹
fastqc -t 5 ./*gz -o ./ #将当前文件夹下所有gz结尾的文件（trim_galore过滤后的文件）进行fastqc处理，生成zip文件和每个文件的html文件，并输出到当前文件夹。
$multiqc ./*zip #对fastqc生成的zip文件整合并生成一个整合后的网页multiqc_report.html。
进入Rstudio: http://172.18.5.194:8889/
直接在Rstudio点击上面生成的两个multiqc_report.html,即可直接查看。
========================================
四：比对
----------------------------------------
下载samtools: 下载方法参考https://www.jianshu.com/p/6b7a442d293f ；使用方法参考https://www.jianshu.com/p/0e3987ee2234
首先，进入github中的samtools链接:https://github.com/samtools/samtools/releases/
找到最下面位于Aseets目录下的安装包samtools.1.9.tar.bz2，然后右键，复制下载地址
wget -c https://github.com/samtools/samtools/releases/download/1.4/samtools-1.4-solo.tar.bz2
tar jxvf samtools-1.4.tar.bz2
cd samtools-1.4.tar.bz2 #进入该目录后发现有一个configure文件
pwd
./configure --prefix=pwd/samtools-1.4 #将samtools的安装路径给prefix
make
make install
./samtools --help #检查是否安装成功
vim ~/.bashrc #添加路径
source ~/.bashrc #更新一下
samtools #检测一下samtools都有哪些命令
总结：源码包的安装一般由3个步骤组成：
配置（configure）
编译（make）
安装 （make install）
安装结束后将自动生成这个包的目录。且所有文档都存在该目录下。使用-prefix选项好处是卸载软件或移植软件。当不需要该软件时只需要简单的删除该安装目录即可把软件卸载干净，移植软件只需拷贝整个目录到另一个机器即可（相同的操作系统）。
直接使用bowtie2比对，需要提前下载索引文件：下载小鼠参考基因组的索引和注释文件，这里常用的mm10。（或者提前下载参考基因组，本次实验中，由于用以下命令下载索引，比对时出现了问题，不知道是不是索引的问题，所以我自己下载了鼠基因组mm39，自己构建索引。由于基因组不大，所以操作起来也很快。）
mkdir reference && cd reference
wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
unzip mm10.zip
以下为测试过程：
cd ~/data/datasource/shengxinpractise/sra/sra && mkdir test_file && touch test1.fq #新建test_file目录和test1.fq文件
cd ~/data/datasource/shengxinpractise/sra/sra && touch test2.fq #新建test2.fq文件
zcat ./SRR2927015_1_val_1.fq.gz |head -10000 > ../test_file/test1.fq #zcat命令查看未解压的压缩包，并将结果输出前10000行到test1.fq文件中（此处选前10000行只是测试需要）
zcat ./SRR2927015_2_val_2.fq.gz |head -10000 > ../test_file/test2.fq #由于是双端测序，故要对另一端的序列进行同样的处理
............


测试无错误，可进行分析了。
mkdir bowtie2_align_sambamba_markdup.sh # 在与clean文件夹并行的目录下创建这个脚本文件
mkdir align #在与clean文件夹并行的文件夹下创建一个align文件夹
vim bowtie2_align_sambamba_markdup.sh
cd ~/data/datasource/shengxinpractise/sra/sra/align
ls ~/data/datasource/shengxinpractise/sra/sra/clean/*_1.fq.gz >1 #将*_1.fq.gz文件（经过trim_galore质控的文件）输出到自动创建的文件夹1中。
ls ~/data/datasource/shengxinpractise/sra/sra/clean/*_2.fq.gz >2 #*_2.fq.gz文件进行同样的处理并输出到文件夹2中。
ls /home/huxl/data/datasource/shengxinpractise/sra/sra/clean/*2.fq.gz | cut -d '/' -f10 | cut -d '_' -f1 >0
#此条命令旨在提取出样本的名称。首先ls输出clean文件夹下的所有*2.fq.gz文件；然后使用cut命令截取字符串，-d '/' 选项指定'/'为分隔符，-f10选项指定在输入文件中即*2.fq.gz文件的所在路径中的第10个域来输出，第10个域为*2.fq.gz，如SRR2927015_2_val_2.fq.gz 文件。）；第二次cut选定分隔符为 '_', -f1选项指定第一个域输出，即输出SRR2927015。
#有关cut 的用法请见https://www.cnblogs.com/dong008259/archive/2011/12/09/2282679.html
#paste 0 1 2 > config.clean #paste命令将0，1,2三个文件夹的相应行用制表符（Tab）连接起来，并输出到config.clean文件中。提示：cat命令也可以进行多个文件内容拼接成一个文件。
sh bowtie2_align_sambamba_markdup.sh #运行这个脚本文件
以上脚本文件运行无误后，再次进入这个脚本文件，开始比对：
vim bowtie2_align_sambamba_markdup.sh
1. 生成一个包含fastq1和fastq2的路径以及样本名称的文件config.clean。
cd ~/data/datasource/shengxinpractise/sra/sra/align
ls ~/data/datasource/shengxinpractise/sra/sra/clean/*_1.fq.gz >1
ls ~/data/datasource/shengxinpractise/sra/sra/clean/*_2.fq.gz >2
ls /home/huxl/data/datasource/shengxinpractise/sra/sra/clean/*2.fq.gz | cut -d '/' -f10 | cut -d '_' -f1 >0
paste 0 1 2 > config.clean

2.进行bowtie2比对
bowtie2_index=~/data/reference/mm10 # 将小鼠参考基因组文件路径赋给变量bowtie2_index
cat config.clean | while read id; #read后加-a选项，则默认只能读出config.clean这个数组的第一列。
# 查看config.clean文件中作为输出|使用read由标准输入读取数据，放入变量id 中，如果读取到的数据非空，就进入循环。
do # 进入循环；参考作者给出的是do echo $id, 我测试过，这样会多输出一次变量id;感觉没必要，就删掉了。直接用do命令即可。
arr=($id) # 将变量id赋值给arr。
fq2=${arr[2]} #将数组arr的第二个列表作为变量赋值给fq2
fq1=${arr[1]} #将arr的第1个列表作为变量赋值给fq1
sample=${arr[0]} #将数组arr的第0个列表作为变量赋值给sample

bowtie2 -p 5 --very-sensitive -X 20000 -x $bowtie2_index -1 $fq1 -2 $fq2 |samtools sort -O bam -@ 5 -o - > ${sample}.raw.bam
#使用bowtie2 命令进行比对；-p 5指定线程数为5；$bowtie2_index输入参考基因组路径； -1 $fq1输入双端测序文件1的路径，-2 $fq2 输入双端测序文件2的路径。samtools -sort对bam文件进行排序，默认是按照序列在fasta文件中的顺序和序列从左往右的位点排序。-O bam指定最终输出文件格式为bam格式，默认是sam格式。-o 后接文件，设置最终排序后的输出文件名。-@ 5指定5个线程。
#实际上，我跑了3个sample,每个sample压缩前文件为4G左右，总共12G，单纯跑这个比对过程，早上十点到晚上十点还没跑完，明天再看看它结束没。

3. sambamba去重。（#samtools去重太麻烦，需要四步。这里选用sambamba去重。）

samtools index ${sample}.raw.bam #为了便于访问bam文件，可以为已经基于坐标排序的bam文件创建索引，-b选项生成以bai索引文件；-c创建csi索引文件.默认参数为-b。（在使用与region参数有关的其他samtool命令时可以使用tabix命令。）

bedtools bamtobed -i ${sample}.raw.bam > ${sample}.raw.bed # #bedtools可以对基因组sam,bed等多种文件格式进行处理，取交集补集等以及格式转换。
#bedtools bambed Usage:bedtools bamtobed [OPTIONS] -i <bam> ；将bam文件转换成bed文件
#此处将 {sample}.raw.bam 文件转换成bed文件并输出到 ${sample}.raw.bed文件中。

samtools flagstat${sample}.raw.bam >${sample}.raw.stat
#统计输入文件的相关数据并将这些数据输出至屏幕显示。每一项统计数据都由两部分组成，分别是QC pass和QC failed，表示通过QC的reads数据量和未通过QC的reads数量。以“PASS + FAILED”格式显示。还可以根据samtools的标志位显示相应的内容，但是这里不做讨论。

sambamba markdup --overflow-list-size 600000 --tmpdir='./' -r ${sample}.raw.bam ${sample}.rmdup.bam
#sambamba markdup命令 默认情况下 marks the duplicates without removing them，可以finding duplicate reads in BAM文件， --overflow-list-size 默认200000，-r选项表示remove instead of mark duplicate. ${sample}.raw.bam 和${sample}.rmdup.bam分别为输入和输出文件。可以是多个输入文件，默认最后一个为输出文件。
#Usage: sambamba-markdup [options] <input.bam> [<input2.bam> [...]] <output.bam>

samtools index ${sample}.rmdup.bam # 为去重之后的输出文件建立索引。

4. calculate %mtDNA;
mtReads=$(samtools idxstats ${sample}.rmdup.bam | grep 'chrM' | cut -f 3)

bowtie2的使用方法：Usage:参考网址：https://blog.csdn.net/herokoking/article/details/77847384
bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r> | --interleaved <i> | -b <bam>} [-S <sam>]
-X 20000 #指定最长插入片段长度为20000；default为500.
-x <bt2-idx> #由bowtie2-build 所生成的索引文件的前缀。首先在当前目录搜寻，然后在环境变量BOWTIE2_INDEXES指定的文件夹中搜寻。
-1 <m1> #双末端测序对应的文件1的路径；
-2 <m2> #双末端测序对应的文件2的路径；
-U <r> #非双末端测序对应的文件，可以为多个文件，用逗号隔开。测序文件的reads长度可以不一样。
-b <bam>} #文件是按read name排序的未对齐BAM文件
-S <sam> #文件按sam文件格式输出；默认是stdout格式

 

未完待续......