以下是根据生信技能树ChIP-seq视频记的笔记，代码来自视频。
1. 安装软件
- 利用以下命令安装所需的软件：
conda install -y sra-tools macs2 deeptools
conda install -y homer meme samtools
conda install -y bowtie bowtie2
conda install -y trim-galore
软件安装在默认目录下 xuyongxin/miniconda3/envs/ChIP-seq/bin
-
安装R包
这里先不安装。 -
查看软件是否安装成功
输入软件名查看help，有帮助命令即安装成功。
2. 下载原始数据
这里使用的数据来自文章：Polycomb associates genome-wide with a specific RNA polymerase II variant, and regulates metabolic genes in ES cells，
网址(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE34520)
- 根据GEO数据库的 accession number GSE34520，在GEO中搜索，点击SRA run selector或者有SRA号，
image.png
当网速较快时，可以找到 accession list，将其下载，用于数据的批量下载。
image.png
- 在主目录下，建立分析所用的文件夹
mkdir -p project/chip-seq
cd project/chip-seq
mkdir {sra,raw,clean,align,peaks,motif,qc}
- 在sra文件夹下建立新文件srr.list，将accession list内容粘贴进去；
- 利用prefetch命令批量下载
如果在conda环境下，对prefetch进行赋值，命令：
prefetch=prefetch
如果不在conda环境下，利用prefetch所在的全路径进行赋值：
prefetch=~/软件所在路径/bin/prefetch
利用srr.list批量下载数据：
cat srr.list | while read id;do (nohup $prefetch $id &);done
这里的$prefetch就是指代上面赋值的“=”后边的prefetch。
利用top可以查看是否正在下载，如果top下没有此任务，可以 用以下命令来搜索任务：
ps -ef | grep pref
默认下载目录为 ~/ncbi/public/sra/，可以更改。
- 因为上述下载老出错，手动安装sratoolkit：
wget https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/current/sratoolkit.current-centos_linux64.tar.gz
tar -xvzf sratoolkit.current-centos_linux64.tar.gz #解压缩
vi ~/.bashrc #添加到环境变量，输入：
export PATH=/data/xuyongxin/biosoft/sratoolkit.2.9.6-1-centos_linux64/bin:$PATH
source ~/.bashrc #激活环境变量
- 使用下面命令批量下载数据：
nohup prefetch --option-file srr.list &
3. 将原始数据转换格式
-
做样本的配置文件，从sra.table提取抗体信息和SRR号信息，包括这两列；
sra.table文件，即下图中的runinfo table文件。
image.png
-
以Tab键为分分割符，我们需要的信息在11，14列，cut出来，然后写perl脚本提取需要的内容。
image.png
cut -f11,14 sra.table > 1_sra_11_14col #冒号前的内容还有
cut -f11,14 sra.table | cut -d ":" -f2 > 1_sra_11_14col #冒号前的内容去掉
Perl脚本如下：
#!/usr/bin/perl -w
my ($infile,$outfile)=@ARGV;
open (IN,"$infile");
open (OUT,">$outfile");
my $header = <IN>;
my %hash;
while (my $line = <IN>){
chomp $line;
my @array = split /\s+/,$line;
my $srr = pop @array;
my $samp = join "_", @array;
if (exists $hash{$samp}){
$hash{$samp} += 1;
}else{
$hash{$samp} = 1;
}
print OUT "$samp\_$hash{$samp}\t$srr\n";
}
close IN;
close OUT;
- 运行以上脚本
perl 1_config_extract.pl 1_sra_11_14col config
- 发现得到的config文件，每行第一个字符都是"_"，去掉：
sed -i 's/^_//' config
至此，得到样本的信息，用于数据格式转换和名称更改。
- 批量转换格式和名称的命令如下，在chip-seq文件夹下运行，config记得拷贝到此处：
analysis_dir=raw
cat config | while read id;
do echo $id
arr=($id)
srr=${arr[1]}
sample=${arr[0]}
nohup $dump -A $sample -O $analysis_dir --gzip --split-3 sra/$srr.sra &
done
注意：上述命令执行是在conda 的ChIP-seq环境下，运行前记得激活环境，且要先对dump进行赋值：
dump=fastq-dump 或者
dump=~/biosoft/sratoolkit.2.9.6-1-centos_linux64/bin/fastq-dump
网友评论