用bowtie2比对blast匹配到lncRNA的序列到基因组上
bowtie2 -p 24 -x ~/RNASEQ/index/bowtie2_index/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.bowtie_index/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.bowtie_index --local -f lncORFs2.txt --al Reads_aligned --un Reads_unaligned | samtools sort -@ 24 -o lncORFs.bam
UCSC上下载GWIPS-viz的bigwig文件
hg19版本的地址是：http://hgdownload.soe.ucsc.edu/gbdb/hg19/bbi/gwipsvizRiboseq.bw
转化为bedgragh文件
./bigwigtobedgragh gwipsvizRiboseq.bw gwipsvizRiboseq.bedGragh
用perl去除低丰度的peak
#!/usr/bin/perl
use strict;
#use warnings;
open (O,$ARGV[0]) or die $!;
open (W,">".$ARGV[1]);
while (<O>) {
chomp;
my @data=split(/\t/,$_);
if ($data[3]>20) {
print W $_."\n";
}
}
bedtools的intersect工具取overlap
bedtools intersect -a lncORFs2.bam -b gwipsvizRiboseq_filted.bed -bed >overlap.bed
overlap注释R code
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
library(clusterProfiler)
library(org.Hs.eg.db)
peaks <- list(gwipsviz="./overlap.bed")
peakAnnoList <- lapply(peaks, annotatePeak, TxDb=txdb, annoDb="org.Hs.eg.db")
anno_tables <- lapply(peakAnnoList, function(i) as.data.frame(as.GRanges(i)))
gwipsviz_anno <- anno_tables$gwipsviz
un_gwipsviz_anno <- gwipsviz_anno[!duplicated(gwipsviz_anno$SYMBOL),]
write.csv(un_gwipsviz_anno,file = "overlap.csv")
