芯片自主注释流程代码

step1-get_fasta
setp2-align
step3-get_gene_position
step4_overlep

参考学习生信菜鸟团和生信技能树的教程
http://www.bio-info-trainee.com/3740.html
http://www.bio-info-trainee.com/3732.html
http://www.bio-info-trainee.com/1469.html

思路

1.获得探针ID和序列

在GEO官网的GPL平台下载gpl
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL21827
这里会看到一个表格，需用getGEO下载GPL文件，从中获取表格中的ID和SEQUENCE，将它构建为fasta格式，第一行为>ID，第二行为sequence。

if(!require(GEOquery)) BiocManager::install("GEOquery")
rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
# 注意查看下载文件的大小，检查数据
f='GPL21827_eSet.Rdata'library(GEOquery)
# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)if(!file.exists('f')){gset <- getGEO('GPL21827', destdir="." )       ## 平台文件save(gset,file='f')   ## 保存到本地
}
load('GPL21827_eSet.Rdata')  ## 载入数据
class(gset)
length(gset)
gset
colnames(Table(gset))
probe2symbol=Table(gset)[,c(1,4)]
all_recs=paste(apply(probe2symbol,1,function(x){paste0('>',x[1],'\n',x[2])}),collapse = '\n') #生成fasta格式的字符串
temp <- tempfile()  ## 编程技巧，把变量写入临时文件~
write(all_recs, temp)

2.比对到参考基因组
首先下载参考基因组和gtf文件，Terminal运行

cd ~/data/project/release1/Genomes/
#Homo_sapiens.GRCh38.dna.toplevel.fa.gz
wget -c ftp://ftp.ensembl.org/pub/release-94/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.toplevel.fa.gz
#Homo_sapiens.GRCh38.94.gtf.gz
wget -c ftp://ftp.ensembl.org/pub/release-94/gtf/homo_sapiens/Homo_sapiens.GRCh38.94.gtf.gz
gunzip *.gz然后用Rsubread构建索引，运行比对library(Rsubread)
# 推荐从ENSEMBL上面下载成套的参考基因组fa及基因注释GTF文件
dir='~/data/project/release1/Genomes/'
ref <- file.path(dir,'Homo_sapiens.GRCh38.dna.toplevel.fa')
buildindex(basename="reference_index",reference=ref)
## 是单端数据
reads <- temp
align(index="reference_index",readfile1=reads,output_file="alignResults.BAM",phredOffset=64)
propmapped("alignResults.BAM")

propmapped函数计算了比对上的探针的百分比。

3.从gtf获取gene位置信息

读取gtf，用getGenePositions和GRanges函数获得一个注释的Grange对象，包括seqid，start，end，strand以及gene_id，明明为my_gr。library(Rsubread)
# 推荐从ENSEMBL上面下载成套的参考基因组fa及基因注释GTF文件
dir='~/data/project/release1/Genomes/'
gtf <- file.path(dir,'Homo_sapiens.GRCh38.82.gtf')
if(!require(refGenome)) install.packages("refGenome")
# create ensemblGenome object for storing Ensembl genomic annotation data
ens <- ensemblGenome()
# read GTF file into ensemblGenome object
read.gtf(ens,useBasedir = F,gtf)class(ens)
# counts all annotations on each seqname
tableSeqids(ens)
# returns all annotations on mitochondria
extractSeqids(ens, 'MT')
# summarise features in GTF file
tableFeatures(ens)
# create table of genes
my_gene <- getGenePositions(ens)
dim(my_gene)# length of genes
gt=my_gene
my_gene_length <- gt$end - gt$start
my_density <- density(my_gene_length)
plot(my_density, main = 'Distribution of gene lengths')library(GenomicRanges)
my_gr <- with(my_gene, GRanges(seqid, IRanges(start, end), strand, id = gene_id))

4.将探针比对到基因组的bam文件坐标信息与来自gtf的坐标信息对应起来

library(Rsamtools)
bamFile="alignResults.BAM"
quickBamFlagSummary(bamFile)
# https://kasperdanielhansen.github.io/genbioconductor/html/Rsamtools.html
bam <- scanBam(bamFile)
bam
names(bam[[1]])
tmp=as.data.frame(do.call(cbind,lapply(bam[[1]], as.character)))
tmp=tmp[tmp$flag!=4,] # 60885 probes
#  intersect() on two GRanges objects.
library(GenomicRanges)
my_seq <- with(tmp, GRanges(as.character(rname), IRanges(as.numeric(pos)-60, as.numeric(pos)+60), as.character(strand), id = as.character(qname)))
gr3 = intersect(my_seq,my_gr)
gr3
o = findOverlaps(my_seq,my_gr)
o
lo=cbind(as.data.frame(my_seq[queryHits(o)]),as.data.frame(my_gr[subjectHits(o)]))
head(lo)
write.table(lo,file = 'GPL21827_probe2ensemb.csv',row.names = F,sep = ',')

从bam文件中读取相同格式的坐标信息，寻找其中有overlap的片段，用cbind建立对应关系。

芯片自主注释流程代码相关推荐

MSM USB插入流程代码分析
点击打开链接代码路径:kernel\msm-3.18\drivers\power\qpnp-smbcharger.c src_detect_handler -->update_usb_stat ...
ATECC508A芯片开发笔记（九）：加密读写508芯片数据的流程及相应设置
目录 ATECC508A芯片开发笔记(九):加密读写508芯片数据的流程及相应设置 1.Encrypted Read 1.1 Standard Encrypted Read Flow 1.2 Simp ...
TI单芯片毫米波雷达1642代码走读（〇）——总纲
前言近年来,自动驾驶行业发展如火如荼,雷达技术也逐渐从军工封闭圈走向了开放的市场. 毫米波雷达具有全天候探测能力,特别是在雨雪雾天气以及夜间都能可靠工作,并且探测距离相对其他车载传感器非常远,对运动 ...
ARM芯片上电启动流程
下图是大多数开发板所有的一个存储单元框架,接下来以此图为基础描述ARM芯片的上电启动流程. 我们首先来了解几个关键词: IROM (Internal ROM):芯片内部固化存储代码的存储器 IRAM ...
完整的芯片反向设计流程原来是这样的！（实例讲解）
完整的芯片反向设计流程原来是这样的!(实例讲解) 作者:时间:2018-02-23来源:网络收藏现代IC产业的市场竞争十分激烈,所有产品都是日新月异,使得各IC设计公司必须不断研发新产品,维持自身企 ...
芯片的设计流程和流片成本
每天都在用,但你知道芯片的设计流程和流片成本吗? 2017-05-10 06:10 来源:半导行业观察芯片,是无数设计工程师们烧死很多脑细胞后产生的作品,完全可以称得上是当代的艺术品.无论是电工们, ...
【嵌入式硬件芯片开发笔记】4-20mA DAC芯片AD5421配置流程
[嵌入式硬件芯片开发笔记]4-20mA DAC芯片AD5421配置流程 16位.串行输入.环路供电.4 mA至20 mA DAC 可用于HART协议相关电路同AD5700配合使用 AD5421的SP ...
【嵌入式硬件芯片开发笔记】HART协议调制解调芯片AD5700配置流程
[嵌入式硬件芯片开发笔记]HART协议调制解调芯片AD5700配置流程 XTAL_EN接地,CLK_CFG的两个引脚由同一个GPIO控制初始时HART_CLK_CFG输出低电平由RTS引脚控制调制 ...
《机器学习实战》萌新入门笔记 ① — — K近邻算法趣味讲解和书本实例详细注释后代码
开坑前言:大一在读新生开启了自己的机器学习之旅,从接触到现在已经有快两个月了,现在回过头为夯实基础,从开始看起<机器学习实战>这本书,即使是最简单的KNN算法,详细看过作者用python实 ...
(十二) 完整注释的代码摘录
title: 完整注释的代码摘录 date: 2019/4/23 20:40:00 toc: true --- 完整注释的代码摘录作者网页 #include <linux/kernel.h&g ...

芯片自主注释流程代码

芯片自主注释流程代码相关推荐

最新文章

热门文章