Lecture 1-10

最近这周帮助组内成员完成一些高重复性基本湿实验、下地收集果实，以及自己的一小部分分子生物学实验外，花了些时间进行了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 的检索和抓取。

NCBI对edirect的介绍

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