数据分析3.1-stringtie 注释新基因-鉴定lncRNA

Posted 2023-05-09

参考技术A Error: input file –e cannot be found!

Aborted (core dumped)

-e选项是只组装参考基因组里有的，不进行新基因的预测，提高组装速度。为什么呢能用呢？

转录组入门：了解参考基因组及基因注释

任务列表

• 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