用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")

GWIPS-viz + hORFeome