这里是佳奥！
让我们从配置环境开始吧！
1 Linux环境
##conda下载安装
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh
##安装
bash Miniconda3-latest-Linux-x86_64.sh
source ~/.bashrc
##配置环境
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
##新建小环境
conda create -n atac python=2 bwa
conda info --envs
source activate atac
##可以用search先进行检索
conda search trim_galore
##保证所有的软件都是安装在atac这个环境下面
conda install -y sra-tools
conda install -y trim-galore samtools
conda install -y deeptools homer meme
conda install -y macs2 bowtie bowtie2
2 R环境
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
install.packages("devtools",
repos="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
library(devtools)
source("https://bioconductor.org/biocLite.R")
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
BiocInstaller::biocLite(c('airway','DESeq2','edgeR','limma'))
BiocInstaller::biocLite(c('ChIPpeakAnno','ChIPseeker'))
BiocInstaller::biocLite('TxDb.Hsapiens.UCSC.hg19.knownGene',
ask=F,suppressUpdates=T)
BiocInstaller::biocLite('TxDb.Hsapiens.UCSC.hg38.knownGene',
ask=F,suppressUpdates=T)
BiocInstaller::biocLite('TxDb.Mmusculus.UCSC.mm10.knownGene',
ask=F,suppressUpdates=T)
##值得注意的是Y叔的包检查会有版本的问题，包括 ChIPseeker
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(ChIPpeakAnno)
library(ChIPseeker)
3 下载sra数据
##测试数据
2-ceLL-1 SRR2927015
2-ce11-2 SRR2927016
2-ce11-5 SRR3545580
2-ce11-4 SRR2927018
touch config.sra
vim config.sra
$ cat config.sra
2-ceLL-1 SRR2927015
2-ce11-2 SRR2927016
2-ce11-5 SRR3545580
2-ce11-4 SRR2927018
cut -d" " -f 2 config.sra > srr.list
$ cat srr.list
SRR2927015
SRR2927016
SRR3545580
SRR2927018
##添加到环境变量
export PATH="$PATH:/home/kaoku/biosoft/sratoolkit/sratoolkit.3.0.0-ubuntu64/bin"
组织好项目
mkdir -p ~/project/atac/
cd ~/project/atac/
mkdir {sra,raw,clean,align,peaks,motif,qc}
cd sra
##把上一级的srr.list移到sra目录
mv ../srr.list ./
批量下载
cat srr.list | while read id
do
nohup prefetch $id &
done
##下载目录在/home/ncbi/sra
QQ截图20221228160006.png
##下载的sra文件
(chipseq) root 11:27:18 /home/kaoku/project/atac/sra
$ ls
SRR2927015.sra SRR2927016.sra SRR2927018.sra SRR3545580.sra srr.list
4 转化为fastq
##fastq-dump
fastq-dump sra/SRR2927015.sra
--gzip :输出文件以gz格式压缩
--split-3 :输入文件为双端测序文件
-A :输出文件名
-O :输出目录
$ cat config.sra
2-ceLL-1 SRR2927015
2-ce11-2 SRR2927016
2-ce11-5 SRR3545580
2-ce11-4 SRR2927018
##循环处理
cat config.sra | while read id;
do echo $id
arr=($id)
srr=${arr[1]}
sample=${arr[0]}
nohup fastq-dump -A $sample -O /home/kaoku/project/atac/raw --gzip --split-3 /home/kaoku/project/atac/sra/$srr.sra &
done
##转化后的fastq文件
(chipseq) root 11:28:23 /home/kaoku/project/atac/raw
$ ls
2-ce11-2_1.fastq.gz 2-ce11-4_1.fastq.gz 2-ce11-5_1.fastq.gz 2-ceLL-1_1.fastq.gz
2-ce11-2_2.fastq.gz 2-ce11-4_2.fastq.gz 2-ce11-5_2.fastq.gz 2-ceLL-1_2.fastq.gz
5 fastq过滤及质控可视化
TrimGalore
##添加TrimGalore到环境变量（依赖cutadapt）
export PATH="$PATH:/home/kaoku/biosoft/trimgalory/TrimGalore-0.6.7"
conda install -c bioconda cutadapt
##trim_galore
-q:设置线程
--phred33 :使用ASCII+33质量值作为Phred得分
--length :去除长度小于参数值的reads
-e :允许的最大误差
--stringency :设置与接头重叠的序列
--paired :对于双端测序文件，正反链质控通过才保留
-o :输出目录
fastq_filel:上一步fastq-dump生成的fastq1文件
fastq_file2:上一步fastq-dump生成的fastq1文件
##构建config
ls *_1.fastq.gz > 1
ls *_2.fastq.gz > 2
paste 1 1 2 > config.raw
##循环过滤脚本
cd /home/kaoku/project/atac/raw
cat config.raw | while read id;
do echo $id
arr=($id)
fq1=${arr[1]}
fq2=${arr[2]}
sample=${arr[0]}
nohup trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired -o /home/kaoku/project/atac/clean $fq1 $fq2 &
done
质控可视化
##raw
cd /home/kaoku/project/atac/qc/raw
fastqc -t 5 /home/kaoku/project/atac/raw/*.gz -o ./
multiqc -n 'raw_fq' ./
##clean
cd /home/kaoku/project/atac/qc/clean
fastqc -t 5 /home/kaoku/project/atac/clean/*.gz -o ./
multiqc -n 'clean_fq' ./
得到了质控过滤后的数据，下一篇我们继续比对。
我们下一篇再见！
网友评论