WES|WGS

巴黎板蓝根

已于 2024-06-24 21:05:17 修改

阅读量43

点赞数

文章标签： linux

于 2023-03-20 13:36:00 首次发布

本文链接：https://blog.csdn.net/JjinYyi/article/details/129665722

版权

这个博主看下测试WGS https://blog.csdn.net/yangl7?type=blog
比较全面的wes流程 https://blog.csdn.net/bio_meimei/article/details/108236406
wes和这个搭配着学习 https://www.jianshu.com/p/cd81bb3be361
WGS处理流程 https://zhuanlan.zhihu.com/p/29485987

0.先检查数据

md5sum -c md5.txt > md5.check
#注1：执行这个命令检查每个文件的md5值是否与文件md5.txt（测序公司给的）中记录的一致。
#注2：极力建议md5.txt中使用相对路径！！！如果路径不对md5值无法匹配检查。
#注3：结果保存在文件md5.check中。对每个文件会给出中文的“正确”，“错误”标识。

在不同批次的文件夹下混合着rna-seq和wes，查找" -T fq.gz"或"-N fq.gz"为RNA-seq(引号内空格处填入星号)

# find -name "*-T*fq.gz" | xargs  -I  '{}'  mv  {}  /public/wgs/JG_PA/RNAseq/FASTQ/
# 教训啊：不要轻易动源文件，要建立软链接，把不同文件夹下符合要求的文件再返回去，没成功，md5对不上，sad
建立软链接：
find dirpath -name "*-T*fq.gz"  -exec realpath {} \; >> ./for_ln_name1.txt
find dirpath -name "*-N*fq.gz"  -exec realpath {} \; >> ./for_ln_name1.txt
# 代码解析："-exec realpath {} \; "用来返回find命令的绝对路径
awk '{print "ln -s",$0,"/public/wgs/JG_PA/RNAseq/FASTQ"}' for_ln_name1.txt >> ln_name.sh

多个批次合并时，可能遇到的问题：

# 1.可能有重复：
# 把文件名输出到不同的txt后，注意文件名要sort成字符串排序的格式
"sort file1.txt" 
"comm -12 file1.txt file2.txt"
# 2.统一文件名：
find ./ -name "*fq*" | while read f; do mv $f ${f/fq/fastq}; done

一、安装GATK+环境 + 构建索引

到这里找对应版本下载 https://github.com/broadinstitute/gatk/releases?page=2

1.下载软件并配置环境

wget https://github.com/broadinstitute/gatk/releases/download/4.1.2.0/gatk-4.1.2.0.zip
unzip gatk-4.1.2.0.zip
cd gatk-4.1.2.0
./gatk
# 配置环境
echo 'export PATH=/home/kelly/biosoft/gatk/gatk-4.1.2.0/:$PATH' >>~/.bashrc
source ~/.bashrc
gatk --help

2.下载参考基因组

“GATK官网推荐使用lftp工具进行访问ftp和下载数据，如果服务器中没有lftp命令，可以提前下载安装，最好使用管理员安装”
官网下载：

location: ftp.broadinstitute.org/bundle
username: gsapubftp-anonymous
password:

或者安装ftp在线下载（上下两步是无法运行的）

sudo apt-get install lftp
#首先lftp，后面跟用户名，然后at符号，ftp服务器地址
lftp ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/
#这里密码是空的，我们直接敲回车即可
mirror hg38 #下载文件夹直接用mirror
#耐心等待
exit #退出当前环境
tree -h #查看当前文件夹内容及大小

因为无sudo权限，手动安装lftp的话也找不到lftp-4.8.4/usr/local文件夹去手动配置（手动配置参考这里https://blog.csdn.net/ct_bao/article/details/103489680，还需要较高版本的gcc，遂放弃。），因此无法ftp连接GATK官网下载； xftp软件也登录不了ftp；GATK官网也进不去（网络不安全报错）

再看这里的时候收到启发 https://www.jianshu.com/p/cd81bb3be361 ：
在这里插入图片描述我们可以截取前半部分下载：

nohup wget -c -r -np -k -L -p ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38
会生成一个“ftp.broadinstitute.org”文件夹，这不行，会断开。
-》结合上下图，一个一个下载这些：
wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G*
wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Axiom*
wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo*
wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp*

在这里插入图片描述

3.bwa构建索引

3.1 安装BWA： https://github.com/lh3/bwa

wget https://github.com/lh3/bwa/releases/tag/v0.7.17
tar -zxvf bwa-v0.7.17.tar.gz
cd bwa-0.7.17/
make
./bwa ## 会弹出帮助信息
vim ~/.bashrc
export PATH=$PATH:/home/liuxin/bwa
source ~/.bashrc

费半天劲还没装成功，还是conda吧

conda create -n WES
conda activate WES
conda install samtools bcftools BWA -y

3.2 构建索引

#先解压human.fa文件，然后构建索引
gzip -d Homo_sapiens_assembly38.fasta.gz ##解压前是849MB，解压后大小是3.1GB
#构建索引，大概3h，格式：bwa index ref.fa -p genome 
##可以不加-p genome，这样建立索引都是以ref.fa为前缀
bwa index Homo_sapiens_assembly38.fasta -p bwa_index/gatk_hg38

BWA命令详解 https://www.jianshu.com/p/19f58a07e6f4
周末要下班回家，就投递任务：

#!/bin/sh
#SBATCH -J genome_index
#SBATCH -p research
#SBATCH -N 1
#SBATCH -n 16 
#SBATCH --error=1.genome_index_err       %指定错误文件输出
bwa index Homo_sapiens_assembly38.fasta -p bwa_index/gatk_hg38
bwa index ref.fa -p genome

将上述代码保存为"1.genome_index.sh"，然后运行：

mkdir bwa_index
sbatch 1.genome_index.sh
scancel 3086269 # 取消作业
scontrol show job 3086269 # 没运行成功，可能因为已完成
sacct -j 3086269 # 查看已结束进程，显示已完成
# 提交后会显示"Submitted batch job 3086269",可以看到"slurm-3086269.out"文件

“1.genome_index_err” 文件打开会显示图1的样子；完成之后会显示图2的样子：