最近这周帮助组内成员完成一些高重复性基本湿实验、下地收集果实,以及自己的一小部分分子生物学实验外,花了些时间进行了Biostar Handbook课程的部分学习,同时,今天跟课题组长沟通了300份核心种质资源的重测序打算,明年准备开展这项工作,以及后面的数据分析。
尽管没有完全看完,但就目前学习到SNP calling这个阶段来说,这本书写的还是非常有实操性,适合于新手拿来进行找感觉。
下面就本书1-10章所涉及的知识点,做一些总结(部分摘抄了生信媛公众号中的内容):
Linux部分
- 常用命令:
#这部分内容都是非?;镜牟僮髅睿ǎ?pwd
mkdir
ls
ls -l
cd
touch <filename> #新建文件
mv
cp
rm -f #强制删除
rm -r #删除文件夹
curl <链接> -o <filename> #用curl命令把URL下载下来,并存储一个叫 filename 的文件中.
curl -O <url> #使得curl 命令去识别url上的文件名(作为下载后的文件名字)
wc -l
head
tail
more
less
grep
grep gene sc.gff | wc -l #管道操作
#cut 在不指定分隔符的时候,默认是tab分割。如果是其他分隔符,可以用-d指定分隔符类型,以逗号分隔(比如csv格式)为例:
cut -d "," -f 1,2,3 features.gff | head
#计算某个词有多少个重复(uniq 加-c后显示重复个数)
cut -f 3 features.gff |sort | uniq -c | head
#模式匹配,是否存在pattern TATA?
gzcat illumina-data.fq.gz | head -10000 | grep 'ATG' --color=always | head
gzcat illumina-data.fq.gz | head -10000 | grep 'TATA' --color=always | head
# 搜索开头为ATG的行
#*在正则表达式里,行开头用^表示
cat SRR519926_1.fastq | egrep "^ATG" --color=always | head
#搜索行末尾是ATG的行
#*在正则表达式里,$表示匹配输入字符串的结束位置
cat SRR519926_1.fastq | egrep "ATG\$" --color=always | head
#*egrep表示扩展正则表达式,用法上grep –E等同于egrep
#*上文中的"\"是转义的作用,这里去除或者保留并不影响匹配结果。
# 搜索TAATA或TATTA模式,这是一个字符范围
cat SRR519926_1.fastq | egrep "TA[A,T]TA" --color=always | head
# 搜索TAAATA 或者 TACCTA, 这些是一些分组
cat SRR519926_1.fastq | egrep "TA(AA|CC)TA" --color=always | head
# 搜索“ TA后有2个或5个A,然后是TA”
cat SRR519926_1.fastq | egrep "TAA{2,5}TA" --color=always | head
# 在任何位置匹配AGATCGG,AGATCGG后可以跟着任何数目的碱基。
cat SRR519926_1.fastq | egrep "AGATCGG.*" --color=always | head
- 环境变量的设置
Linux中的PATH是一个字符串变量,里面记录了一系列目录,当运行一个程序时,Linux在这些目录下进行搜寻编译链接。
对PATH进行编辑时,其格式为:
PATH=$PATH:<PATH1>:<PATH2>:<PATH3>:------:<PATHN>
export命令用于设置或显示环境变量,可新增,修改或删除环境变量,供后续执行的程序使用。export的效力仅及于该次登陆操作。
语法:
export [-fnp][变量名称]=[变量设置值]
查看PATH可使用:
echo $PATH
在Linux中进行环境变量的添加可用三种方法(以添加blast+2.7.1为例):
方法一:
直接用export命令:
export PATH=$PATH:~/src/ncbi-blast-2.7.1+/bin
配置完后可以通过echo $PATH查看配置结果。
生效方法:立即生效
有效期限:临时改变,只能在当前的终端窗口中有效,当前窗口关闭后就会恢复原有的path配置
用户局限:仅对当前用户
方法二:
通过修改.bashrc文件:
vim ~/.bashrc
#在最后一行添上:
export PATH=$PATH: ~/src/ncbi-blast-2.7.1+/bin
#或者在之前的PATH变量声明后添加:
export PATH=$PATH:~/src/edirect:~/src/ncbi-blast-2.7.1+/bin
生效方法:(有以下两种)
1、关闭当前终端窗口,重新打开一个新终端窗口就能生效
2、输入“source ~/.bashrc”命令,立即生效
有效期限:永久有效
用户局限:仅对当前用户
方法三:
通过修改profile文件:
vim /etc/profile
#在最后一行添上:
export PATH=$PATH: ~/src/ncbi-blast-2.7.1+/bin
#或者在之前的PATH变量声明后添加:
export PATH=$PATH:~/src/edirect:~/src/ncbi-blast-2.7.1+/bin
生效方法:(有以下两种)
1、关闭当前终端窗口,重新打开一个新终端窗口就能生效
2、输入“source /etc/profile”命令,立即生效
有效期限:永久有效
用户局限:对所有用户
一个栗子:
# 将~/bin 文件夹加到PATH
echo 'export PATH=~/bin:$PATH' >> ~/.bashrc
source ~/.bashrc
# 在~/bin生成fastqc快捷方式
ln -s ~/src/FastQC/fastqc ~/bin/fastqc
软件使用部分
1. Entrez Direct
edirect 全称 EntrezDirect,是 NCBI 开发的用于 linux 命令行界面的快速检索和下载工具,使用它可以很方便的在服务器上进行 Entrez 的检索和抓取。
1.1 下载与安装:
wget ftp://ftp.ncbi.nlm.nih.gov/entrez/entrezdirect/edirect.tar.gz
tar zxvf edirect.tar.gz
sh ~/app/edirect/setup.sh
echo "export PATH=\$PATH:\$HOME/src/edirect" >> $HOME/.bashrc
1.2 主要工具与用途
工具名称 | 用途 |
efetch | 下载 NCBI 数据库中的记录和报告并以相应格式打印输出 |
einfo | 获取目标结果在数据库中的信息 |
elink | 对目标结果在其他数据库中比配结果 |
epost | 上传 UIDs 或者 序列登记号 |
esearch | 在 Entrez 中执行搜索命令 |
efilter | 对之前的检索结果进行过滤或限制 |
xtract | converts XML into a table of data values. |
nquire | sends a URL request to a web page or CGI service. |
运行实例:
#抓取nucleotides数据。
esearch -db nucleotide -query PRJNA257197 | efetch -format fasta > ebola.fasta
# 以GenBank格式获取数据
esearch -db nucleotide -query PRJNA257197 | efetch -format gb > ebola.gb
# 通过locus号来获取一条序列。
efetch -db nucleotide -id KM233090 -format fasta > KM233090.fa
#一次性获取多个id下的序列。
efetch -db nucleotide -id KM233090,KM233066,KM233113.1 -format fasta > multi.fa
#从数据中只获取一个子集。
efetch -db nucleotide -id KM233090,KM233066,KM233113.1 -format fasta -seq_start 1 -seq_stop 10 > multi.fa
2. SRA toolkit
2.1 下载与安装:
curl -O https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.8.2/sratoolkit.2.8.2-ubuntu64.tar.gz
tar xzvf sratoolkit.2.8.2-ubuntu64.tar.gz
echo export PATH=$PATH:~/src/edirect:~/src/sratoolkit.2.8.2-ubuntu64/bin >> ~/.bashrc
source ~/.bashrc
2.2 使用
#获取文件
prefetch SRR1553610
#该文件放置在~/ncbi文件夹下
#用程序fastq-dump来把文件拆包(针对双端测序pair-end)
fastq-dump --split-files SRR1553610
# 如何下载多个文件?创建一个含有SRR runs的文件。
echo SRR1553608 > sra.ids
echo SRR1553605 >> sra.ids
prefetch --option-file sra.ids
#只拆包sra.ids里指定的文件,用到sed (字符流编辑器) 工具来提取文件并替换其中的模式“SRR”
cat sra.ids | sed 's/SRR/fastq-dump --split-files SRR/' | bash
#sed中的"s"表示substitute,用于替换紧跟/ /内的内容
3. FastQC与trimmomatic
3.1下载与安装
#安装FastQC:
cd ~/src
curl -O http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.zip
unzip fastqc_v0.11.5.zip
cd FastQC
chmod +x fastqc #设置可执行。
# 安装trimmomatic
cd ~/src
curl -O http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.32.zip
unzip Trimmomatic-0.32.zip
cd Trimmomatic-0.32
java -jar trimmomatic-0.32.jar #运行需要java
# 我们写一个脚本来替代我们运行。
echo '#!/bin/bash' > ~/bin/trimmomatic
echo 'java -jar ~/src/Trimmomatic-0.32/trimmomatic-0.32.jar $@' >> ~/bin/trimmomatic
chmod +x ~/bin/trimmomatic
# 看,成功了!
trimmomatic
3.2 运行
#!/bin/bash
# 这个脚本会把目录下的所有fastq文件进行修剪。
for name in *.fastq
do
echo "*** currently processing file $name, saving to trimmed-$name"
trimmomatic SE -phred33 $name trimmed-$name SLIDINGWINDOW:4:30 MINLEN:35 TRAILING:3
echo "*** done"
done
# 现在,运行fastqc报告。
# 你可以在这添加更多的命令。
fastqc trimmed-*.fastq
# 查看FastQC 和Trimmomatic的污染序列和接头序列文件。
ls ~/src/FastQC/Configuration/
# 文件内容涵盖一些已知的接头序列。
more ~/src/FastQC/Configuration/adapter_list.txt
# 文件内容涵盖已知的污染序列。
more ~/src/FastQC/Configuration/contaminant_list.txt
# 进行接头切除时,我们需要找到含有接头序列的文件。
# 你可以自己创建自己的接头文件,或者使用Trimmomatic自带的
# 我们来创建一个Illummina的接头文件吧。
echo ">adapter" > adapter.fa
echo "AGATCGGAAGAG" >> adapter.fa
# 我们来进行质量和接头的修剪。
trimmomatic SE SRR519926_1.fastq trimmed.fq ILLUMINACLIP:adapter.fa:2:30:10 SLIDINGWINDOW:4:30 TRAILING:30 ILLUMINACLIP:adapter.fa:2:30:5
# Trimmomatic自带了一些接头,我们可以来用一下。
# 依据文件内容不同,可能有不同的操作模式。
# 这是一个palindromic接头。 Trimmomatic通过接头名字来确定操作模式。
#*留意你的安装版本号。如果安装了更新版本,需要对应改一下路径名称。
ln -s ~/src/Trimmomatic-0.32/adapters/TruSeq3-PE.fa
# 这是一个扩展的palindromic接头文件,它也列出了单个接头。
ln -s ~/src/Trimmomatic-0.32/adapters/TruSeq3-PE-2.fa
# 这些是一些经典的接头。
ln -s ~/src/Trimmomatic-0.32/adapters/TruSeq3-SE.fa
trimmomatic SE SRR519926_1.fastq trimmed.fq ILLUMINACLIP:adapter.fa:2:30:10 SLIDINGWINDOW:4:30 TRAILING:30 ILLUMINACLIP:TruSeq3-SE.fa:2:30:5
# 在双端(paired end)模式下运行时,两个测序文件需要被同时过滤和剪切。
# 命令行混乱随之而来。
trimmomatic PE SRR519926_1.fastq SRR519926_2.fastq trim_1P.fq trim_1U.fq trim_2P.fq trim_2U.fq ILLUMINACLIP:TruSeq3-PE-2.fa:2:10:10 SLIDINGWINDOW:4:30 TRAILING:20
# 一个稍微简单的选择是,加一个 -baseout