Draw GO bubble Picture with results of dataframe format of clusterProfiler (根据clusterProfiler的dataframe结果,绘制气泡图)

library(ggplot2)
options(stringsAsFactors=F)
dz <- read.table("test.txt",sep="\t",header=T,check.names=F)
head(dz)

n <- 1
for (i in (dz$GeneRatio)){
    j=strsplit(i,'\\/')
    dd <- unlist(j)
    value <- (as.numeric(dd[1])/as.numeric(dd[2]))
    dz$GeneRatio[n] <- sprintf("%.2f",value)
    n=n+1
    }


pp <- ggplot(dz,aes(GeneRatio,Description))
pp <- pp + geom_point()
pbubble <- pp + geom_point(aes(size=Count,color=qvalue))
tiff(file="GO_dotplot_new.tiff",width=25,height=30,units="cm",compression="lzw",bg="white",res=600)
pbubble + scale_colour_gradient(low="green",high="red") + labs(color=expression(qvalue),size="Count",x="GeneRatio",y="GO results",title="Results of GO enrichment analysis") + facet_grid(ONTOLOGY~.,scale="free")
dev.off()

发表在 R | 留下评论

deal with results of clusterProfiler and draw the personal picture(根据需要处理clusterProfiler的结果)

library(devtools)
install_github("YuLab-SMU/clusterProfiler.dplyr")
library(clusterProfiler.dplyr)
ddd <- filter(kk_all_readable,Count>=4) #根据自己需要对数据过滤
tiff(file="KEGG_dotplot.tiff",width=25,height=30,units="cm",compression="lzw",bg="white",res=600)
dotplot(ddd, x="GeneRatio",color = "pvalue")
dev.off()

参考链接:clusterProfiler.dplyr

发表在 R | 留下评论

TCGA_ATAC-seq 比对文件整理 (Adjusting alignment file of ATAC-seq in TCGA)

#!/usr/bin/perl -w
use strict;
use warnings;
use Cwd;
use File::Copy;

my $usage=<<USAGE;
Usage:
    perl $0 filename
    将该目录下的所有文件放在一个叫filename的文件夹下
USAGE
if(@ARGV==0){die $usage};

my $file1=$ARGV[0];
my $dir=getcwd;
my @filename;

opendir(DIR,$dir) || die $!;
my @temp=readdir(DIR);
for my $i (0..$#temp){
    if ((-d "$dir/$temp[$i]") & ($temp[$i] ne $file1)){
        opendir(DIR2,"$dir/$temp[$i]") || die $!;
        my @temp1=readdir(DIR2);
        for my $j (0..$#temp1){
            if((-f "$dir/$temp[$i]/$temp1[$j]") & (($temp1[$j])=~/.bam$/) | ($temp1[$j]=~/.bai$/) ){
                push(@filename,"$dir/$temp[$i]/$temp1[$j]");
            }
        }
    }
}
close(DIR2);
close(DIR);

#print @filename;

my $lujing="$dir/$file1";

if(!(-e $lujing)){
    mkdir ($lujing) || die "can not create $file1 doctory";
}

for my $i (0..$#filename){
    move("$filename[$i]","$lujing");
}
发表在 Perl, TCGA | 留下评论

use wget when prefetch disabled (prefetch无法下载时,使用wget下载)

less -SN nohup.out | tail -n 40 | grep https:// | perl -e 'while(<>){chomp($_);$_=~/(https:.*)/;print $1."\n"}'| sed 's/^/wget\t-c\t-t\t0\t/g' | perl -e 'while(<>){chomp;my @arr=split(/\t/,$_);my @temp=split(/\//,$arr[4]);my $id=pop(@temp);$id=~/(\S+)\.\d/;print $_."\t"."-O ".$1.".sra"."\n";}' | sed 's/\t/ /g' >> download_wget.sh
cat download_wget.sh | while read id ;do (nohup $id &);done
发表在 Linux, NGS_analysis, Perl | 留下评论

Deal with documents by Perl

#!/usr/bin/perl -w
use strict;
use warnings;

my $usage=<<USAGE;
Usage:
    perl $0 filename_with_clinical.txt /path/to/file
USAGE
if(@ARGV==0){die $usage};

my $file1=$ARGV[0];
my $file2=$ARGV[1];
my @filename;

opendir(DIR,$file2) || die $!;
@filename=readdir(DIR);
splice(@filename,0,2);
close(DIR);

my $tcgaid;

open(RF,$file1) || die $!;
while(my $line=<RF>){
    chomp($line);
    my @arr=split(/\t/,$line);
    if(($. % 2) == 0){
        $tcgaid="$arr[2]-T2-$arr[3]";
    }else{
        $tcgaid="$arr[2]-T1-$arr[3]";
    }
    for my $i (0..$#filename){
        if($arr[0] eq $filename[$i]){
            rename ("$file2/$filename[$i]","$file2/$tcgaid.bw");
        }
    }
}
close(RF);
发表在 Perl | 留下评论

handle many SAM files (批量处理SAM文件)

#sam2bam,bam2sort,bam2index
screen -S yourname
for i in `ls *.sam`;
do
j=${i/\.sam/};
echo "samtools view -b -S $j.sam > $j.bam
samtools sort $j.bam sorted_$j
samtools index sorted_$j.bam" >> $j.sh
done

for i in `ls *.sh`;do echo "bash $i" >> total.sh;done

ParaFly -c total.sh -CPU 9
发表在 Linux | 留下评论