探针重注释到lncRNA

先看曾老师完全用R骚操作的教程，反正我的小本本是带不起来这样的操作的：

• 搜狗微信搜索，第一篇就是

首先去affy官网下载相应的fasta文件：
hg-u133-plus

然后用bowtie2比对到hg38基因组：

source activate wes2
ref=/home/qmcui/database/reference/index/bowtie/hg38
bowtie2 -x $ref -f HG-U133_Plus_2.probe_fasta -S hgu133plus2.sam
samtools view -b -S hgu133plus2.sam > hgu133plus2.bam
bamToBed -i hgu133plus2.bam > hgu133plus2.bed

完成后的bed文件：
前3列是基因组坐标，第4列是包含探针id的信息，第6列是正负链

去gencode下载lncRNA的gtf文件：

lncRNA gtf

wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_29/gencode.v29.long_noncoding_RNAs.gtf.gz
gunzip gencode.v29.long_noncoding_RNAs.gtf.gz

然后用bedtools通过取交集的方式对前面比对好的探针进行注释：

intersectBed -a hgu133plus2.bed -b gtf2bed/gencode.v29.annotation.gtf -wa -wb >hgu133plus2.txt

红色框内是一行信息，黄色框内色是hgu133plus2.bed文件信息，其余是gencode.v29.annotation.gtf文件信息

接下来需要对信息进行筛选：

cat hgu133plus2.txt | awk '$9=="gene"&&$6==$13{print $4,$16,$18, $20}'> hgu133plus2lnc.txt

只要第9列是'gene'，且第6列和第13列标明的正负链方向一致的行。

image.png
image.png

把文件下载到本地，接下来就进入R中操作了。

options(stringsAsFactors = F)
library(data.table)
library(stringr)
rm(list=ls())

probe2lnc = fread('./hgu133plus2lnc.txt')
probe2lnc = probe2lnc[, -5]
colnames(probe2lnc) = c('probeID', 'EnsID', 'type', 'Symbol')
#对探针id的字符串进行提取
probe2lnc$probeID = word(probe2lnc$probeID, 3, sep = ':')

tmp = unique(probe2lnc)#去除重复行
#对probeID进行计数，如果出现两次或两次以上
#说明一个探针映射到多个基因，要舍去该probe
tmp$ProbeNum = table(tmp$probeID)[tmp$probeID]
hgu133plus2lnc = tmp[tmp$ProbeNum == 1,-5]

save(hgu133plus2lnc, file = 'hgu133plus2lnc.rdata')

得到的结果：

image.png

---

生信技能树公益视频合辑：学习顺序是linux，r，软件安装，geo，小技巧，ngs组学！
B站链接：https://m.bilibili.com/space/338686099
YouTube链接：https://m.youtube.com/channel/UC67sImqK7V8tSWHMG8azIVA/playlists
生信工程师入门最佳指南：https://mp.weixin.qq.com/s/vaX4ttaLIa19MefD86WfUA
学徒培养：https://mp.weixin.qq.com/s/3jw3_PgZXYd7FomxEMxFmw

其实我还是喜欢简书的：
生信人的Linux20题