以下是根据生信技能树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号,
当网速较快时,可以找到 accession list,将其下载,用于数据的批量下载。
- 在主目录下,建立分析所用的文件夹
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文件。
-
以Tab键为分分割符,我们需要的信息在11,14列,cut出来,然后写perl脚本提取需要的内容。
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