刘小泽写于2018年，发送于2019.6.13
相信你已经磨拳擦掌，想要自己分析一波原始的人类基因组数据，并自己挑出变异基因，这感觉很像一个侦探，逐步深入案情，十分吸引人！
1 计算资源准备
最低配置Linux 8线程+16G+1T硬盘【因为原始数据很大，参考的各个数据库也很大】
查看linux 下计算机配置：
CPU：cat /proc/cpuinfo | grep "process" | wc -l
内存：free -g
硬盘：df -h
2 原始数据下载(放在~/project/kpgp_wgs/raw_data下)
选取韩国人Korean Personal Genome Project中的一个人的测序数据KPGP-00001
nohup wget -c -r -nd -np -k -L ftp://ftp.kobic.re.kr/pub/KPGP/2015_release_candidate/WGS/KPGP-00001
关于wget这几个参数：-c 是断点续传；-r是递归下载，因为这里提供的下载地址相当于一个文件夹，下面有很多数据；-nd （no-directories ）不创建目录；-np（no-parent） 不要追溯到上一级目录；-k（ –convert-links） 转换非相对链接为相对链接【相当于linux中的绝对路径，能够直接访问下载的；因为网页中为了让链接容易记住，设置了大量的相对连接，而我们下载相对链接是访问不进去的】；-L（-L, –relative ）仅仅跟踪相对链接
这些参数设置也不是一成不变的，一般需要设置断点-c、递归-r、-nd
3 参考数据下载
所有的下载之前都要新建好明确的文件夹，比如用户是在home目录下的，一开始空空如也。新建几个目录mkdir biosoft(存放软件)、mkdir project（存放项目）、mkdir reference（存放参考数据）
3.1 参考基因组（放在refernce中）
cd ~/reference
# human genome 19（hg19）
mkdir -p genome/hg19 && cd genome/hg19
for i in $(seq 1 22) X Y M; #指定染色体号下载
do echo $i;
# 一般下载UCSC的数据
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz;
done
gunzip *.gz
for i in $(seq 1 22) X Y M; # 指定染色体号拼成整体
do cat chr${i}.fa >> hg19.fasta
done
rm -fr chr*.fa
# hg38一样的道理
# 参考基因组大小压缩包是900M左右，解压完是3G
小疑点，大问题～基因组版本
人类基因组主要存放在三个数据库NCBI、Ensembl、UCSC
NCBI：ftp://ftp.ncbi.nlm.nih.gov/genomes/Homo_sapiens/ARCHIVE/
Ensembl：ftp://ftp.ensembl.org/pub/【然后找release版本，最新是release-93】
UCSC：http://hgdownload.cse.ucsc.edu/goldenPath/hg19
【UCSC的hg系列，是目前使用频率最高的基因组。最新是hg38】
3.2 下载好的基因组需要构建索引（放在reference中）
练习使用bowtie2,hisat2和bwa这3个主流比对软件
软件安装先一律使用
conda install -p /YOUR/PATH/ soft_name -y
-p制定安装路径；-y允许安装
3.2.1 bowtie2
cd ~/reference
mkdir -p index/bowtie2 && cd index/bowtie2
nohup time bowtie2-build --threads 20 ~/reference/genome/hg19/hg19.fa hg19 1>hg19.bowtie2.log 2>&1 &
#time计算运行时间；--threads设置线程
# 1>hg19.bowtie2.log 新建hg19.bowtie2.log，将运行信息存入；
# 2>&1 2代表错误信息，若有报错，将错误信息覆盖到1的文件中
#hg38参考基因组同理
3.2.2 hisat2
cd ~/reference
mkdir -p index/hisat2 && cd index/hisat2
nohup time hisat2-build -p 20 ~/reference/genome/hg19/hg19.fa hg19 1>hg19.hisat2.log 2>&1 &
#-p设置线程
#hg38参考基因组同理
#也可以下载现有的：ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/
#下载后的文件夹中还包括make_hg19.sh这个程序，运行可以直接下载hg19基因组
3.2.3 bwa
cd ~/reference
mkdir -p index/bwa && cd index/bwa
nohup time bwa index -a bwtsw -p hg19 ~/reference/genome/hg19/hg19.fa 1>hg19.bwa_index.log 2>&1 &
#-a指定算法bwtsw；-p指定输出名前缀
#hg38参考基因组同理
3.3 下载基因组注释文件
下载之前先建好存放数据文件夹，先下载gencode数据库中的注释文件。hg38修改了hg19的一些错误，但是在提取指定区域时容易出现一个基因对应两个“分身”位置【见参考4】
这里就有一个问题：GTF和GFF一样吗？
GTF是在GFF的基础上发展来的，二者都是
\t分隔的9列文件；GFF包含的信息更多，可以有染色体、基因、转录本信息，而GTF(Gene transfer format)主要描述基因和转录本
mkdir -p ~/reference/gtf/gencode && cd ~/reference/gtf/gencode
## GRCh38 https://www.gencodegenes.org/releases/current.html
mkdir GRCh38_hg38 && cd GRCh38_hg38
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gtf.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.2wayconspseudos.gtf.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.long_noncoding_RNAs.gtf.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.polyAs.gtf.gz
## GRCh37
mkdir ~/reference/gtf/gencode/GRCh37_hg19 && cd ~/reference/gtf/gencode/GRCh37_hg19
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/GRCh37_mapping/gencode.v28lift37.annotation.gtf.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/GRCh37_mapping/gencode.v28lift37.metadata.HGNC.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/GRCh37_mapping/gencode.v28lift37.metadata.EntrezGene.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/GRCh37_mapping/gencode.v28lift37.metadata.RefSeq.gz
然后按照我们的需要提取出一些指定区域，比如蛋白质编码区，长链非编码RNA区、假基因区(pseudogene是与正常基因相似，但丧失正常功能的DNA序列)以及hg19、hg38各自的全基因位点区域信息
## 提取GRCh37(hg19)全基因信息
cd ~/reference/gtf/gencode/GRCh37_hg19
zcat gencode.v28lift37.annotation.gtf.gz | perl -alne '{next unless $F[2] eq "gene" ;/gene_name \"(.*?)\";/; print "$F[0]\t$F[3]\t$F[4]\t$1"}'> allgene.hg19.position
对这一行代码的解读：
zcat是在不解压缩文件的前提下查看文件，结果通过管道符|传给后面命令；perl接受了gtf文件中一行一行的信息，要统计各个基因所在染色体、起始位置、终止位置、基因名这四种信息，就要先对gtf文件做解读
处理tab分割的文件的每一行内容，
perl -alne是必备佳品：
例如：while(<>){chomp;@F=split /s+/,$_;print "$F[2]\n"}也就相当于
perl -alne 'print $F[2]～多么简洁！
-e打开单行输入模式；-n打开循环模式；-l对输入内容自动chomp，对输出内容自动添加换行符；-a自动分隔模式
那么这里的{next unless $F[2] eq "gene意思就是如果第二列（Feature信息）不等于“gene”就看下一行；
/gene_name \"(.*?)\";/是将gene_name后面引号中(注意引号要用转义一下)的部分存入待引用变量中，那么存入多少呢？这里的(.*?)就起作用了：.是除了\n的任意字符，*是取之前字符的0个或者n个。perl默认使用贪婪模式匹配，也就是如果不限制，他能一直匹配到结尾，因此?就是限制贪婪模式的，这三个连用，就命令perl，只准匹配到后一个引号，不能继续了哈！
接下来就是打印输出了，print "$F[0]\t$F[3]\t$F[4]\t$1"将第一列、第四列、第五列再加上前面待引用部分【用$1调取】打印出来，中间用分隔符分隔
要知道提取什么位置，就必须要先了解gtf的内容，除了之前说的那几列，有一个很重要的知识点：gene type的类型重复多次出现，有什么规律吗？
##统计hg19编码蛋白的基因信息
zcat gencode.v28lift37.annotation.gtf.gz | grep protein_coding | perl -alne '{next unless $F[2] eq "gene"; /gene_name \"(.*?)\";/;print "$F[0]\t$F[3]\t$F[4]\t$1"}'>hg19_protein_coding.position
##统计hg38的长链非编码RNA信息
zcat gencode.v28.long_noncoding_RNAs.gtf.gz| perl -alne '{next unless $F[2] eq "gene"; /gene_name \"(.*?)\";/;print "$F[0]\t$F[3]\t$F[4]\t$1"}'>hg38_lncRNA.position
##统计hg38的假DNA信息(第三列type信息全是transcript)
zcat gencode.v28.2wayconspseudos.gtf.gz | perl -alne '{next unless $F[2] eq "transcript"; /gene_name \"(.*?)\";/;print "$F[0]\t$F[3]\t$F[4]\t$1"}'> hg38_pseudos.position
##统计hg38的全基因信息
zcat gencode.v28.annotation.gtf.gz | perl -alne '{next unless $F[2] eq "gene"; /gene_name \"(.*?)\";/;print "$F[0]\t$F[3]\t$F[4]\t$1"}'>allgene.hg38.position
结果得到了hg19的编码蛋白基因信息
要查找基因信息只需要grep基因名就好，标准基因名可在https://www.genecards.org/中查找，例如要查找抑癌基因TP53和乳腺癌基因BRCA1的信息，就可以
grep TP53 hg19_protein_coding.position
grep BRCA1 hg19_protein_coding.position
欢迎关注我们的公众号～_～
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com
Welcome to our bioinfoplanet!
网友评论