一个简单进行SNP分析的实战例子（优秀学徒成果展）
准备工作
1
修改提示符
具体详见：shell界面颜值知多少 (颜值即正义)
1
保存到环境变量
mkdir biosoft package
mkdir ~/bisoft/mybin
vim .bashrc
export PATH=$PATH:/home/yjzhang/biosoft/mybin
安装软件
首先参考：生物信息学常见1000个软件的安装代码！ (初学者先这样摸索)
1
安装fastq-dump (sratoolkit套件)
cd package
wget -c http://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.9.1-1/sratoolkit.2.9.1-1-ubuntu64.tar.gz
tar -zxvf sratoolkit.2.9.1-1-ubuntu64.tar.gz -C ~/biosoft/
cd ~/biosoft/sratoolkit.2.9.1-1-ubuntu64/bin
ln -s ~/biosoft/sratoolkit.2.9.1-1-ubuntu64/bin/fastq-dump ~/biosoft/mybin/
2
安装Entrez Direct软件套件
curl ftp://ftp.ncbi.nlm.nih.gov/entrez/entrezdirect/edirect.zip -O
unzip edirect.zip -d ~/biosoft/
## 不能执行,应该没有安装上，找教程发现还有要执行的步骤
./setup.sh
#添加到环境变量
echo "export PATH=\${PATH}:/home/yjzhang/biosoft/edirect" >> $HOME/.bashrc
安装后有帮助信息出现，应该是安装成功了，但下载数据失败，由于时间紧迫，目前没找原因，
查询conda的软件库有这个软件，直接转战conda安装
conda 软件查询网址
https://bioconda.github.io/recipes.html
3
安装miniconda
wget -c https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
chmod +x Miniconda3-latest-Linux-x86_64.sh
./Miniconda3-latest-Linux-x86_64.sh
##更改镜像配置
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 snp
source activate snp #启动新环境
source deactivate #关闭环境
4
conda install entrez-direct
5
conda install readseq
6
conda install bwa
7
conda install samtools
8
conda install freebayes
具体就不罗列了，总之是严格按照生信技能树的教学视频来的，自己去B站观看：
下载原始数据
从埃博拉项目中获取多个数据集
项目编号：PRJNA257197
这个时候需要学习生信技能树的：解读SRA数据库规律一文就够
1
第一种方法（下载的完整数据，比较大）
NCBI搜索，下载该项目的SraRunInfo.csv文件，FileZilla导入到服务器
mkdir sra
cd sra
cat SraRunInfo.csv | cut -f 10 -d ',' | grep SRR | head -10 > runinfo.csv
cat runinfo.csv | xargs -i echo wget -c {} \$ >>sra.sh
bash sra.sh
##下载完后进行文件转换
cd ~/snp/sra
cat SraRunInfo.csv | cut -f 1 -d ',' | grep SRR | head -10 > sraids.txt
#批量转换
vim fastq.sh
#! /bin/bash
#
bin=fastq-dump
dir='/home/yjzhang/snp/raw/'
cat $1 |while read id
do
arr=($id)
echo "Processin sample $arr"
$bin --gzip --split-3 -O $dir /home/yjzhang/snp/sra/$arr
done
bash fastq.sh sraids.txt
2
第二种下载方法
##下载完后进行文件转换
下载的数据不完整，每个文件只下载了20000条数据，方便练习，节省时间
mkdir ~/snp/sra/fastq
cat SraRunInfo.csv | cut -f 1 -d ',' | grep SRR | head -10 > sraids.txt
cat sraids.txt | awk ' { print "fastq-dump -X 20000 --split-files -O ~/snp/sra/fastq " $1 } ' > get-data.sh
bash get-data.sh
制作参考数据
cd ~/snp/ref
efetch -id KM034562 -format gb -db nucleotide > KM034562.gb
readseq -format=FASTA -o ~/snp/ref/KM034562.fa KM034562.gb
readseq -format=GFF -o ~/snp/ref/KM034562.gff KM034562.gb
比对
1
构建索引
bwa index ~/snp/ref/KM034562.fa
samtools faidx /KM034562.fa
2
比对
mkdir ~/snp/align
cd ~/snp/align
vim bwa.sh #新建脚本
NAME=$1
READ1=~/snp/sra/fastq/"$NAME"_1.fastq
READ2=~/snp/sra/fastq/"$NAME"_2.fastq
REFS=~/snp/ref/KM034562.fa
GROUP="@RG\tID:$NAME\tSM:$NAME\tLB:$NAME"
BAM=~/snp/align/$NAME.bam
bwa mem -R $GROUP $REFS $READ1 $READ2 2> logfile.txt | samtools view -b - | samtools sort -o $BAM -
echo $?
samtools index $BAM
echo "*** Created alignment file $BAM"
samtools flagstat $BAM | grep proper
cat ~/snp/sra/sraids.txt | awk ' { print "bash bwa.sh " $1 } ' > align-bwa.sh
bash align-bwa.sh
找SNP
cd ~/snp/align
freebayes -p 1 -f ~/snp/ref/KM034562.fa ~/snp/align/*.bam > freebayes.vcf
把bam文件，vcf文件，和参考序列导入IGV查看
随便找一个基因
samtools view -h SRR1972880.bam KM034562:9547-9637 |samtools sort -o 9547.bam -
这个流程就是生信技能树的WES学徒课程，学完马上就能follow啦，还等什么呢，赶快试试看