任务列表

  • 1.在UCSC下载hg19参考基因组;
  • 2.从gencode数据库下载基因注释文件,并且用IGV去查看感兴趣的基因的结构,比如TP53,KRAS,EGFR等等。
  • 3.截图几个基因的IGV可视化结构
  • 4.下载ENSEMBL,NCBI的gtf,也导入IGV看看,截图基因结构
  • 5.了解IGV常识

在UCSC下载hg19参考基因组

hg19、GRCH38、 ensembl75这3种基因组版本应该是大家见得比较多的了,国际通用的人类参考基因组,其实他们储存的是同样的fasta序列,只是分别对应着三种国际生物信息学数据库资源收集存储单位,即NCBI,UCSC及ENSEMBL各自发布的基因组信息而已。有一些参考基因组比较小众,存储的序列也不一样,比如BGI做的炎黄基因组,还有DNA双螺旋结构提出者沃森(Watson)的基因组,还有2016年发表在nature上面的号称最完善的韩国人做的基因组。前期我们先不考虑这些小众基因组,主要就下载hg19和hg38,都是UCSC提供的,虽然hg38相比hg19来说,做了很多改进,优点也不少,但因为目前为止很多注释信息都是针对于hg19的坐标系统来的,我们就都下载了,正好自己探究一下。也顺便下载一个小鼠的最新版参考基因组吧,反正比对也就是睡个觉的功夫,顺便分析一下结果,看看比对率是不是很低。

mkdir rna_seq/data/reference && cd rna_seq/data/reference
mkdir -p genome/hg19 && cd genome/hg19
# nohup wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz &
# nohup 是永久执行,& 是指在后台运行。nohup COMMAND & 这样就能使命令永久的在后台执行
nohup axel http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz &
tar zvfx chromFa.tar.gz
cat *.fa > hg19.fa
rm chr*.fa

从gencode数据库下载基因注释文件,并且用IGV去查看感兴趣的基因的结构

下载基因注释文件

官网:http://www.gencodegenes.org/releases/26lift37.html

wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_26/GRCh37_mapping/gencode.v26lift37.annotation.gtf.gz
gzip -d gencode.v26lift37.annotation.gtf.gz

 

下载安装IGV、BEDtool

官网:http://software.broadinstitute.org/software/igv/download(下载 Binary Distribution 版本)

wget https://github.com/arq5x/bedtools2/releases/download/v2.26.0/bedtools-2.26.0.tar.gz
tar -zxvf bedtools-2.26.0.tar.gz
cd bedtools2
make

 

截图几个基因的IGV可视化结构

批量截图:TP53,KRAS,EGFR

grep -w 'gene' gencode.v26lift37.annotation.gtf | grep -w 'TP53' | cut -f 1,4,5 >> gene.bed
grep -w 'gene' gencode.v26lift37.annotation.gtf | grep -w 'KRAS' | cut -f 1,4,5 >> gene.bed
grep -w 'gene' gencode.v26lift37.annotation.gtf | grep -w 'EGFR' | cut -f 1,4,5 >> gene.bed
~/biosoft/bedtools2/bin/bedtools igv -i gene.bed > Bach_sanpshot.txt

grep是一个多用途的文本搜索工具,linux中使用非常频繁,并且使用很灵活,可以是变量,也可以是字符串。最基本的用法有以下两种:

  • 1.搜索内容中无空格,可以直接执行grep命令,比如:grep pass a.txt,表示在a.txt文件中搜索pass所在的行
  • 2.如果搜索内容中有空格,则需要使用单引号或者双引号把搜索内容引起来,比如:grep "hello all" a.txt或者grep 'hello all' a.txt,如果不加单双引号,则提示错误,无法识别,因为不加引号,直接grep hello all a.txt,表示在all和a.txt中搜索hello,这肯定是不对的

grep -w option file:精确搜索,可以说准确性搜索,比如:grep -w b* a.txt:此命令执行时,*不会默认为任何字符,只表示字面意思,就是一个*字符

管道命令操作符:”|”,它仅能处理经由前面一个指令传出的正确输出信息,也就是 standard output 的信息,对于 stdandard error 信息没有直接处理能力。然后,传递给下一个命令,作为标准的输入 standard input

cut 命令从文件的每一行剪切字节、字符和字段并将这些字节、字符和字段写至标准输出。如果不指定 File 参数,cut 命令将读取标准输入。必须指定 -b、-c 或 -f 标志之一。使用 -f 选项提取指定字段

下载ENSEMBL,NCBI的gtf

axel ftp://ftp.ensembl.org/pub/grch37/release-89/gtf/homo_sapiens/Homo_sapiens.GRCh37.87.gtf.gz
axel ftp://ftp.ensembl.org/pub/grch37/release-89/gtf/homo_sapiens/Homo_sapiens.GRCh37.87.chr.gtf.gz
axel  ftp://ftp.ncbi.nlm.nih.gov/genomes/Homo_sapiens/ARCHIVE/ANNOTATION_RELEASE.105/GFF/ref_GRCh37.p13_top_level.gff3.gz
axel ftp://ftp.ncbi.nlm.nih.gov/genomes/Homo_sapiens/ARCHIVE/ANNOTATION_RELEASE.105/GFF/ref_GRCh37.p13_scaffolds.gff3.gz