查看原文
其他

一个简单进行SNP分析的实战例子(优秀学徒成果展)

生信技能树 生信技能树 2022-06-06

准备工作

 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啦,还等什么呢,赶快试试看


独家福利


如果需要组装自己的服务器;代办生物信息学服务器

如果需要帮忙下载海外数据(GEO/TCGA/GTEx等等),点我?

如果需要线下辅导及培训,看招学徒

如果需要个人电脑:个人计算机推荐

如果需要置办生物信息学书籍,看:生信人必备书单

如果需要实习岗位:实习职位发布

如果需要售后:点我


您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存