BOL: All

Create a heatmap with R

Rahul Nayak — Mon, 31 Jul 2017 08:45:58 -0500

bio <- read.csv("ppg2008.csv", sep=",")

bio <- bio[order(bio$PTS),]
row.names(bio) <- bio$Name
bio <- bio[,2:20]

bio_matrix <- data.matrix(bio)
bio_heatmap <- heatmap(bio_matrix, Rowv=NA, Colv=NA, col = brewer.pal(9, "Blues"), scale="column", margins=c(5,10))


##
#Sample DATA
#Name  ,G,MIN,PTS,FGM,FGA,FGP,FTM,FTA,FTP,3PM,3PA,3PP,ORB,DRB,TRB,AST,STL,BLK,TO,PF
#Genome1 ,79,38.6,30.2,10.8,22,0.491,7.5,9.8,0.765,1.1,3.5,0.317,1.1,3.9,5,7.5,2.2,1.3,3.4,2.3
#Genome2 ,81,37.7,28.4,9.7,19.9,0.489,7.3,9.4,0.78,1.6,4.7,0.344,1.3,6.3,7.6,7.2,1.7,1.1,3,1.7
#Genome3,82,36.2,26.8,9.8,20.9,0.467,5.9,6.9,0.856,1.4,4.1,0.351,1.1,4.1,5.2,4.9,1.5,0.5,2.6,2.3


library(ggplot2)
bio <- read.csv("seeTNF_Final", sep="\t")
row.names(bio) <- bio$Contig
bio <- bio[,2:256]
data=as.matrix(bio)
head(data)
#Rcolorbrewer palette
library(RColorBrewer)
coul = colorRampPalette(brewer.pal(8, "PiYG"))(25)
#heatmap(data)

# Use 'scale' to normalize (right)
heatmap(data, scale="column")
#heatmap(data, scale="column", col = coul)

Construct a heatmap with perl

Neel — Fri, 07 Jul 2017 04:35:13 -0500

#!/bin/env perl
use strict;
use warnings;

use Imager::Heatmap;

my $hmap = Imager::Heatmap->new(
    xsize  => 640,        # Image width
    ysize  => 480,        # Image height
    xsigma => 10,         # Sigma value of X-direction
    ysigma => 10,         # Sigma value of Y-direction
);

# @point_datas should be: ( [ x1, y1, weight1 ], [ x2, y2, weight2 ] ... )
my @point_datas = ( [ 10, 20, 50 ], [ 20, 40, 70 ] );

# Add point datas to construct density matrix
$hmap->insert_datas(@point_datas); 

# After adding datas, get heatmap as Imager instance.
my $img = $hmap->draw;

# create png file in current directory
$img->write( file => './hm.png' );

Format the fasta header

Jit — Sun, 18 Jun 2017 14:12:39 -0500

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

    my %id2seq = ();
    my $id = 0;
    open F,"$ARGV[0]",or die $!;
    while(){
        chomp;
        if($_ =~ /^>(.+)/){
		$id++;
		s/\[//g; s/\]//g;
		my @val = split /\s+/, $_;
		my @newval=@val;
		shift @newval;
		print ">$id [ gene=$val[0] ] [ protein=@newval ]\n";
            #$id = $1;
        }else{
            print "$_\n";
        }
    }
close F;

Compress and decompress the sequence with perl

Jit — Fri, 16 Jun 2017 09:39:26 -0500

use strict;
use warnings;

my @char;
while () {
@char = split //;
}

comp(\@char);
#---------------------
my $com= "r0a3m4a4j0";
my @com = split //, $com;
dcomp (\@com);

#dcomp sub here
sub dcomp {
my ($com_ref)=@_;
my @com=@$com_ref;
my $car;
for (my $aa=0; $aa<=$#com; $aa++) {
	if ($com[$aa]!~ /\D/) { print "$car" x ($com[$aa]);} else {print "$com[$aa]"; $car = "$com[$aa]";}
}
print "\n";
}

#comp sub here
sub comp {
my ($char_ref)=@_;
my @char=@$char_ref;
my $cnt="";
for (my $aa=0; $aa<=$#char; $aa++) {
	if ($char[$aa+1]) {
		if ($char[$aa] eq $char[$aa+1]) { $cnt++; }  
		elsif ($cnt) { print "$char[$aa]"."$cnt"; undef $cnt;}
		else {print "$char[$aa]";}
	}
}
print "\n";
}
__DATA__
raaaammmmmaaaaaj

Extract fasta sequence from a multifasta file with fasta header Ids

Abhimanyu Singh — Tue, 13 Jun 2017 07:49:37 -0500

#!/usr/bin/perl

use strict;
use warnings;

#Usage: perl   

my $list = shift @ARGV;
my $fasta = shift @ARGV;
my $out = shift @ARGV;
my %select;

open LIST, "$list" or die;
while () {
    chomp;
    s/>//g;
    $select{$_} = 1;
}
close LIST;

$/ = "\n>";
open OUT, ">$out" or die;
open FASTA, "$fasta" or die;
while () {
    s/>//g;
    my ($id) = split (/\n/, $_);
    print OUT ">$_" if (defined $select{$id});
}
close FASTA;
close OUT;

Download the gff files from NCBI using bash script/command

Rahul Nayak — Thu, 08 Jun 2017 08:17:11 -0500

#!/bin/bash

# Download the genome from NCBI using command

# Create a Directory
mkdir genome_gff
cd genome_gff

# Look for genome assembly summary and extract the URL
# USER need to provide the right summary file to curl  
# Commentline if you are not interested in that genome set
# -for fungi
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/fungi/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.gff.gz|' > genomic_file_fungi

# -for bacteria
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/genbank/bacteria/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCA_.+)|\1\2/\2_genomic.gff.gz|' > genomic_file_bacteria

# -for plant 
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/plant/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.gff.gz|' > genomic_file_plant 

# -for archaea
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/archaea/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.gff.gz|' > genomic_file_archaea

# -for protozoa
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/protozoa/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.gff.gz|' > genomic_file_protozoa

# -for vertebrate_mammalian
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_mammalian/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.gff.gz|' > genomic_file_vertebrate_mammalian

# -for vertebrate_other
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_other/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.gff.gz|' > genomic_file_vertebrate_other

# -for invertebrate
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/invertebrate/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.gff.gz|' > genomic_file_invertebrate

# -for viral
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/viral/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.gff.gz|' > genomic_file_viral

#Read the uerl from file and download

FILES=$(pwd)/*
for f in $FILES
do
  echo "Processing $f file..."
  filename=$(basename "$f")
  extension="${filename##*.}"
  filename="${filename%.*}"
  # Create a directory with appending G
  mkdir "GFF$filename"
  cd "GFF$filename"
  # take action on each file. $f store current file name
  head -n 4 $f > $f.head
  wget --input $f.head
  gunzip *.gz
  #cat $f
  cd ..
done

Read a tab delimited file and search with perl

Jit — Tue, 06 Jun 2017 09:18:05 -0500

use strict;
use warnings;
use Data::Dumper;
use Text::CSV;
use IO::Handle;

my $file = "/home/urbe/Tools/Alienomics_v0.1/Alienomics/output/intermediate_files/rRNA/refGene.megablast";
open my $fh, "<", $file or die "$file: $!";
my $csv = Text::CSV->new ({
      sep_char=>"\t",
      binary    => 1, # Allow special character. Always set this
      auto_diag => 1, # Report irregularities immediately
      });

my $fh = IO::Handle->new_from_fd( $fh, 'r' );

while ( not $fh->eof ) {
  my $row = $csv->getline( $fh );
print "$row->[0]\n";
  warn Dumper $row; # To see the structure
}

Unzip all the genome file and remove all fasta header except first one

Jit — Tue, 30 May 2017 10:34:25 -0500

#!/bin/bash

gzip -d *.gz

FILES=$(pwd)/*
for f in $FILES
do
  echo "Processing $f file..."
  	if [[ $f =~ \.fna$ ]];
	then awk ' /^>/ && FNR > 1 {next} {print $0} ' $f | sed '/^>/{s/ /_/g}' > $f.fa 
  	#then sed '1!{/^\>/d;}' $f > $f.fa
  	else echo "this file is not right file"
  	fi
  	#cat $f
done

Download the genome from NCBI using bash script/command

Jit — Tue, 30 May 2017 08:11:33 -0500

#!/bin/bash

# Download the genome from NCBI using command

# Create a Directory
mkdir genome
cd genome

# Look for genome assembly summary and extract the URL
# USER need to provide the right summary file to curl  
# Commentline if you are not interested in that genome set
# -for fungi
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/fungi/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.fna.gz|' > genomic_file_fungi

# -for bacteria
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/genbank/bacteria/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.fna.gz|' > genomic_file_bacteria

# -for plant 
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/plant/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.fna.gz|' > genomic_file_plant 

# -for archaea
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/archaea/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.fna.gz|' > genomic_file_archaea

# -for protozoa
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/protozoa/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.fna.gz|' > genomic_file_protozoa

# -for vertebrate_mammalian
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_mammalian/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.fna.gz|' > genomic_file_vertebrate_mammalian

# -for vertebrate_other
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_other/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.fna.gz|' > genomic_file_vertebrate_other

# -for invertebrate
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/invertebrate/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.fna.gz|' > genomic_file_invertebrate

# -for viral
curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/viral/assembly_summary.txt' | awk '{FS="\t"} !/^#/ {print $20} ' | sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/.+/)(GCF_.+)|\1\2/\2_genomic.fna.gz|' > genomic_file_viral

#Read the url from file and download

FILES=$(pwd)/*
for f in $FILES
do
  echo "Processing $f file..."
  filename=$(basename "$f")
  extension="${filename##*.}"
  filename="${filename%.*}"
  # Create a directory with appending G
  mkdir "G$filename"
  cd "G$filename"
  # take action on each file. $f store current file name
  wget --input $f
  #cat $f
  cd ..
done

#Reference
#ftp://ftp.ncbi.nlm.nih.gov/pub/factsheets/HowTo_Downloading_Genomic_Data.pdf

Extract fasta sequence from a multifasta file with coordinates

Jit — Tue, 30 May 2017 05:35:07 -0500

#!/usr/bin/perl

use Bio::DB::Fasta;

#USAGE
perl extractFASTAwithSIZE.pl finalSample_filtered.fa 0 1000 > aaaaaa.fa

my $fastaFile = shift;
my $querySizeST = shift;
my $querySizeED = shift;

my $db = Bio::DB::Fasta->new( $fastaFile );
my @ids = $db->get_all_primary_ids;

foreach my $id (@ids) {
    my $sequence = $db->seq($id, $querySizeST => $querySizeED);
    if  (!defined( $sequence )) {
            die "Sequence $seq not found. \n" 
    }   
    print ">$id\n", "$sequence\n";
}

__END__

use Bio::DB::Fasta;

  # Create database from a directory of Fasta files
  my $db       = Bio::DB::Fasta->new('/path/to/fasta/files/');
  my @ids      = $db->get_all_primary_ids;

  # Simple access
  my $seqstr   = $db->seq('CHROMOSOME_I', 4_000_000 => 4_100_000);
  my $revseq   = $db->seq('CHROMOSOME_I', 4_100_000 => 4_000_000);
  my $length   = $db->length('CHROMOSOME_I');
  my $header   = $db->header('CHROMOSOME_I');
  my $alphabet = $db->alphabet('CHROMOSOME_I');

  # Access to sequence objects. See Bio::PrimarySeqI.
  my $seq     = $db->get_Seq_by_id('CHROMOSOME_I');
  my $seqstr  = $seq->seq;
  my $subseq  = $seq->subseq(4_000_000 => 4_100_000);
  my $trunc   = $seq->trunc(4_000_000 => 4_100_000);
  my $length  = $seq->length;

  # Loop through sequence objects
  my $stream  = $db->get_PrimarySeq_stream;
  while (my $seq = $stream->next_seq) {
    # Bio::PrimarySeqI stuff
  }

  # Filehandle access
  my $fh = Bio::DB::Fasta->newFh('/path/to/fasta/files/');
  while (my $seq = <$fh>) {
    # Bio::PrimarySeqI stuff
  }

  # Tied hash access
  tie %sequences,'Bio::DB::Fasta','/path/to/fasta/files/';
  print $sequences{'CHROMOSOME_I:1,20000'};