Recherche de tag: R


Parseur de csv basique [Python]

04.07.2018     qcavaille      python3 parser csv 

  Parseur de csv basique :
- Délimiteur réglable
- Pas de limite de colonnes, mais nécessite des lignes "complètes" ET une ligne d'en-têtes de colonnes
- Retourne un générateur, théoriquement plus léger en mémoire. Mais itérable une seule fois !
import csv as L
def readCSV(file):
    with open(file) as csvfile:
        readCSV = L.reader(csvfile, delimiter=";")
        firstRow=True
        columnName=[]
        for row in readCSV:
            if firstRow:
                firstRow=False
                for col in row:
                    columnName.append(col)
            else:
                aRow={}
                i=0
                for col in columnName:
                    aRow[col]=row[i]
                    i+=1
                yield aRow
2/5 - [3 ratings]




reverse complement a sequence (probably the best solution or near of them) [C++]

04.07.2018     Natir      rev_comp singleton 

  A powerful reverse complement algorithm, encapsulating it in a singleton to avoid unnecessary multiple allocation.

#include <string>

class rev_comp
{
public:
    static std::string run(std::string seq)
    {
	rev_comp::build_instance();
	
	auto first = seq.begin(), last = seq.end();
        
	while(true)
        {
		if(first == last || first == --last)
        	{
        	    if(seq.length() % 2)
        		    *first = rev_comp::complement[*first];
        	    return seq;
        	}
        	else
        	{
        	    *first = rev_comp::complement[*first];
        	    *last = rev_comp::complement[*last];
        	    std::iter_swap(first, last);
        	    ++first;
            }
        }
    }

protected:
    static rev_comp* _instance;
    static void build_instance()
    {
        if(_instance == nullptr)
            _instance = new rev_comp();
    }
    
    rev_comp()
    {
        this->complement['A'] = 'T';
        this->complement['T'] = 'A';
        this->complement['C'] = 'G';
        this->complement['G'] = 'C';
        this->complement['a'] = 't';
        this->complement['t'] = 'a';
        this->complement['c'] = 'g';
        this->complement['g'] = 'c';
    }
        
    ~rev_comp(){}
    
    
private:
    static char complement['t' + 1];
};

char rev_comp::complement['t' + 1];
rev_comp* rev_comp::_instance = nullptr;
5/5 - [2 ratings]




Diminution de la taille d'une matrice [R]

04.07.2018     gdevailly      matrice interpolation heatmap 

  Ploter une grosse matrice sous forme d'heatmap peut prendre du temps.
Il peux être judicieux de réduire la taille de la matrice avant de la ploter, en moyennant des groupes de cellules adjacentes (ou en prenant le maximum, lorsque la matrice est 'sparse').
C'est notamment utile lorsqu'il y a beaucoup plus de cellules que de pixels sur les écran actuels (par exemple, ploter 80 000 régions génomiques sur un écran HD de 1080 pixels de hauts).

Utilisation:

> bigmat <- matrix(1:64, ncol = 8)
> bigmat
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 9 17 25 33 41 49 57
[2,] 2 10 18 26 34 42 50 58
[3,] 3 11 19 27 35 43 51 59
[4,] 4 12 20 28 36 44 52 60
[5,] 5 13 21 29 37 45 53 61
[6,] 6 14 22 30 38 46 54 62
[7,] 7 15 23 31 39 47 55 63
[8,] 8 16 24 32 40 48 56 64
> # default usage
> redim_matrix(bigmat, target_height = 4, target_width = 3)
[,1] [,2] [,3]
[1,] 10.0 30.0 50.0
[2,] 11.5 31.5 51.5
[3,] 13.0 33.0 53.0
[4,] 15.0 35.0 55.0
> # changing aggregating function
> redim_matrix(bigmat, target_height = 4, target_width = 3, summary_func = function(x) max(x, na.rm = TRUE))
[,1] [,2] [,3]
[1,] 19 43 59
[2,] 20 44 60
[3,] 22 46 62
[4,] 24 48 64
> # multicore
> redim_matrix(bigmat, target_height = 4, target_width = 3, n_core = 2)
[,1] [,2] [,3]
[1,] 10.0 30.0 50.0
[2,] 11.5 31.5 51.5
[3,] 13.0 33.0 53.0
[4,] 15.0 35.0 55.0
# reduce matrix size, using a summarizing function (default, mean)
redim_matrix <- function(
    mat,
    target_height = 100,
    target_width = 100,
    summary_func = function(x) mean(x, na.rm = TRUE),
    n_core = 1
    ) {

    if(target_height > nrow(mat) | target_width > ncol(mat)) {
        stop("Input matrix must be bigger than target width and height.")
    }

    seq_height <- round(seq(1, nrow(mat), length.out = target_height + 1))
    seq_width  <- round(seq(1, ncol(mat), length.out = target_width  + 1))

    # complicate way to write a double for loop
    do.call(rbind, parallel::mclapply(seq_len(target_height), function(i) { # i is row
        vapply(seq_len(target_width), function(j) { # j is column
            summary_func(
                mat[
                    seq(seq_height[i], seq_height[i + 1]),
                    seq(seq_width[j] , seq_width[j + 1] )
                    ]
            )
        }, 0.0)
    }, mc.cores = n_core))
}
5/5 - [1 rating]




Réarrangement de l'ordre des branches d'un dendrogramme [R]

04.07.2018     gdevailly      hclust dendrogram optimal leaf ordering OLO seriation clustering heatmap 

  Lors d'un clustering hierarchique, l'agencement des branches à chaque nœud est laissé libre (telle branche peut être agencé à gauche ou à droite).
Il peut en résulter des visualisations de heatmaps clusterisée parfois pas très belle.
Heureusement, il est possible de réarranger les branches d'un clustering hiérarchique pour aboutir à de jolies images, tel qu'illustré dans ce billet de blog:
http://nicolas.kruchten.com/content/2018/02/seriation/

Cette petite fonction remplacera donc avantageusement la fonction hclust() dans vos heatmaps clusterisées.
# perform hclust, then aply optimal leaf ordering of the dendrogram branches
require(seriation) # you should install.packages("seriation") if error message

hclust_olo <- function(mdist, ...) {
    myClust <- hclust(mdist, ...)
    myOlo <- seriation::seriate(mdist, method = "OLO")
    seriation::permute(myClust, myOlo)
}
5/5 - [1 rating]




Parseur de fichier gff3 basique [Python]

08.07.2018     propan2one      gff3 parser python3 

  La fonction gff3parser() permet de parser un fichier au format .gff3 et de retourner un dictionnaire au format :

{ numéro de l'ORF : [start, stop, strand ], ... }

Les print en commentaire permettent d'afficher le contenue du gff3 sur le terminal

contentGFF = ["seqid", "source", "type", "start", "end", "score", "strand", "phase", "attributes"]
def gff3parser (filename, dictOrf):
	 """ parser de gff3 """
	with open(filename) as f:
		for line in f:
			if line.startswith("#"):
		                continue
		content = line.strip().split("\t")
		if len(content) == len(contentGFF):
			value = [content[3], content[4], content[6]]
			#print(value)
			#print(content[8])
			numORF = content[8].strip().split(";")
			#print(numORF)
			for product in numORF:
				if re.search(r"^product=ORF\d{0,3}", product):
					keys = product[11:]
					dictOrf[keys] = value

		return {}

# Usage :
orf = {}
gff3parser(gffname,orf) #avec gffname le nom du fichier gff3 à parser
1/5 - [1 rating]




Python cProfile as a decorator [Python]

12.07.2018     Plopp      python cprofile profiling profiler decorator 

  Quickly profile your functions with cProfile by using it as a decorator
import cProfile

def profileit(func):
    def profiled_func(*args, **kwargs):
        profile = cProfile.Profile()
        print("Profiling " + "\033[95m" + func.__name__ + "\033[0m") # colors
        try:
            profile.enable()
            result = func(*args, **kwargs)
            profile.disable()
            return result
        finally:
            profile.print_stats(sort="time")
    return profiled_func



##########
# Example #
##########

@profileit
def times2():
   for i in range(500000):
       j = i * 2
       j = add_one(j)

def add_one(n):
    return n+1

times2()
5/5 - [1 rating]




reverse complement a sequence [Perl]

17.07.2018     jbrunier      revcomp perl 

  a simple reverse complement function
my %C = (  'N' => 'N', 'A' => 'T', 'T' => 'A', 'U' => 'A', 'C' => 'G', 'G' => 'C',
                   'M' => 'K', 'K' => 'M', 'R' => 'Y', 'Y' => 'R', 'W' => 'W', 'S' => 'S',
                   'V' => 'B', 'B' => 'V', 'H' => 'D', 'D' => 'H', 'X' => 'X',
                   'n' => 'n', 'a' => 't', 't' => 'a', 'u' => 'a', 'c' => 'g', 'g' => 'c',
                  'm' => 'k', 'k' => 'm', 'r' => 'y', 'y' => 'r', 'w' => 'w', 's' => 's',
                  'v' => 'b', 'b' => 'v', 'h' => 'd', 'd' => 'h', 'x' => 'x', '-' => '-' );


=head2 RevCompl

 Title   : RevCompl
 Usage   : RevCompl($sequence)
 Function: Reverse complement a sequence
 Example : 
 Returns : the reverse complement sequence
 Args    : $sequence	the sequence to reverse complement
 Notes   : 
 
=cut
sub RevCompl {
    join( '', map { $C{$_} } reverse( split( '', shift @_ ) ) );
}
0/5 - [0 rating]




Visualisations d'ACP interactives et colorées ! [R]

31.07.2018     gdevailly      r PCA ACP plotly interactive plots 

  Visualiser des ACP (PCA in english) est une tâche assez commune, mais qui peut-être laborieuse si vous êtes pointilleux comme moi ou ZaZo0o (cf https://bioinfo-fr.net/representer-facilement-une-acp-avec-r-et-ggplot2).

Je vous présente ici une méthode pour faire des ACP interactives (avec zoom, masquage des groupes, noms des points en mouseover, etc.).

Il vous faut installer dplyr, cowplot et plotly depuis CRAN.
library(cowplot)
library(plotly)
library(dplyr)

plotPCA <- function(
    my_pca,
    colorby = NULL,
    xaxs = "PC1", yaxs = "PC2",
    title = NULL,
    palette = function(x) rainbow(x, s = 0.6),
    continuous_color = FALSE,
    ...
) {
    
    # obtaining % of variance explained by each axis
    varxaxs <- (summary(my_pca)$importance["Proportion of Variance", xaxs] * 100) %>% round(digits = 1)
    varyaxs <- (summary(my_pca)$importance["Proportion of Variance", yaxs] * 100) %>% round(digits = 1)
    
    myplot <- tibble(
        myxaxs = my_pca$x[, xaxs],
        myyaxs = my_pca$x[, yaxs],
        texts = rownames(my_pca$x),
        colors = if (is.null(colorby)) {NA} else {colorby}
    ) %>%
        ggplot(aes(x = myxaxs, y = myyaxs, color = colors, text = texts)) +
        geom_point() +
        labs(x = paste0(xaxs, " (", varxaxs, "%)"), y = paste0(yaxs, " (", varyaxs, "%)"), title = title, color = NULL) +
        background_grid()
    
    if (continuous_color) {
        myplot <- myplot + scale_color_gradientn(colors = palette(64))
    } else {
        myplot <- myplot + scale_color_manual(values = palette(n_distinct(colorby)))
    }
    
    ggplotly(myplot, tooltip = "text", ...)
    
}

# generating data for the example
metadata <- tibble(
    group = rep(c("group1", "group2", "group3"), each = 50),
    names = paste("sample", 1:150),
    latent_variable = c(
        rnorm(100) + 4,
        rnorm(50)
    )
)
my_data <- matrix(rnorm(150 * 15), nrow = 150)
my_data[1:100, 1:10] <- my_data[1:50, 1:10] + 3 # PC1 axis differences
my_data[1:50, 11:15] <- my_data[1:50, 11:15] + 2 # PC2 differences 
rownames(my_data) <- metadata$names

# basic example (click on legend to hide/display specific groups):
plotPCA(prcomp(my_data), colorby = metadata$group)

# changing axis and title
plotPCA(prcomp(my_data), colorby = metadata$group, xaxs = "PC2", yaxs = "PC3", title = "PC1 not shown")

# changing color palette:
plotPCA(prcomp(my_data), colorby = metadata$group, palette = viridis::viridis)

# continuous colors:
plotPCA(prcomp(my_data), colorby = metadata$latent_variable, palette = viridis::viridis, continuous_color = TRUE)

# saving as html:
p <- plotPCA(prcomp(my_data), colorby = metadata$group)
htmlwidgets::saveWidget(p, "my_great_pca.html", selfcontained = TRUE)


0/5 - [0 rating]




Compter le nombre de lignes en fonction du contenu d'une colonne [Bash]

07.08.2018     anicolas      comptage occurence tsv 

  
Exemple
Vous voulez compter le nombre d'occurence de chaque type, à partir du fichier tabulé suivant :
id value type
1 5698 typeA
2 569 typeB
3 3658 typeC
4 532 typeB
5 123 typeA

Vous obtenez :
typeA 2
typeB 2
typeC 1

Usage :
Le fichier tabulé doit avoir une ligne avec les noms de colonnes
Pour le lancer
./count_file.sh fichier_tabulé.txt numero_de_colonne_à_compter
#!/bin/bash
#s'il manque un argument, le script n'est pas lance
if [ $# -ne 2 ]; then echo 'Syntaxe: count_file.sh DE_file_with_annot.txt column_number_to_count'; exit 66; fi

#ecrit un fichier temporaire sans les noms de colonnes
tail -n +2 $1 >> $1.temp
#recupere le prefixe du nom de fichier pour ecrire le nom du fichier de sortie
prefix=$(basename $1 .txt)
output=$prefix"_column_num_"$2"_count.txt"
#creation du fichier de sortie et du fichier temporaire de sortie
temp="temp.txt"
touch $output
touch $temp

#lecture du fichier tabule ligne par ligne
while read line
do
#recupere la valeur de la colonne a compter
current_value=`echo "$line" | awk -F $"\t" '{print $'$2'}'`
#cherche dans le fichier de sortie si cette valeur existe. Si elle existe, augmente d'un le compteur pour cette valeur. Si elle n'existe pas, elle est cree et son compteur est initie a 1. Ecrit l'ancien fichier de sortie + les nouvelles infos dans un fichier temporaire
awk -F $"\t" -v value="$current_value" 'BEGIN{ OFS = FS ; found = "no"} {if( $1 == value ) {print $1, $2+1 ; found = "ok"} else {print $0}} END {if (found == "no") {print value, 1} }' $output >> $temp
#remplace le fichier de sortie par le fichier temporaire
mv $temp $output
done<$1.temp

#efface le fichier tabule sans les noms de colonnes
rm $1.temp
1/5 - [1 rating]




Test de Mac-Nemar par regroupement pour analyse diachronique par mailles [R]

27.08.2018     Remy_Moine      MacNemar boucle treemap 

  Permet de calculer pour différents regroupements d'individus statistiques (ici des mailles de présence-absence) les valeurs du test de Mac Nemar.
Demande en entrée un tableau avec le nom du regroupement, le nombre d'individus statistiques où le caractère observé est absent, apparait, disparait et se maintient entre t et t+1.

Packages additionels: "treemap", "RColorBrewer"

R version 3.4.4
OS: Linux Xubuntu 16.04
setwd()

library(treemap)
library(RColorBrewer)

tabl #contient vos données  en entrée classe data.frame

## Préparation des vecteurs nécessaires au fonctionnement de du script:

Somme<-rep(0,nrow(tabl)) # stock le nombre total d'individus statistiques
stat_mcnemar<-rep(0,nrow(tabl)) # stock la VTZ² pour chaque regroupement
pvalue_mcnemar<-rep(0,nrow(tabl)) # stock la p.value pour chaque regroupement
Applicable<-rep(NA,nrow(tabl)) # stock la condition d'appliation du test
Linkimg<-rep(0,nrow(tabl)) # stock le lien vers la sortie graphique

## Boucle de calcul des paramètres du test et des sorties graphiques par regroupement

for (i in c(1:nrow(tabl))){

  # Vérifie que le test de Mac-Nemar est applicable
  Somme[i]<-sum(tabl[i,c(1:4)])
  Applicable[i]<- tabl[i,3]+tabl[i,2]

  # Produit une matrice pour le test de Mac-Nemar
  mat<-matrix(c(tabl[i,4],tabl[i,2],tabl[i,3],tabl[i,1]),ncol=2)

  # Stockage des résultats du test avec signalement des regroupements où il n'est pas applicable
  if (Applicable[i] < 20) stat_mcnemar[i]<- "test non applicable"
  else stat_mcnemar[i]<-mcnemar.test(mat)$statistic

  if (Applicable[i] < 20) pvalue_mcnemar[i]<- "test non applicable"
  else pvalue_mcnemar[i]<-mcnemar.test(mat)$p.value

  # Export du graphique résumant la répartition des mailles sur le site et stockage du chemin d'accès dans le tableau
  sign<-NULL
  if (mcnemar.test(mat)$p.value == "NaN") sign<- c("NA") else
    if (mcnemar.test(mat)$p.value < 0.01) sign<-c("***") else
      if (mcnemar.test(mat)$p.value < 0.05) sign<-c("**") else
        if (mcnemar.test(mat)$p.value < 0.1) sign<-c("*") else sign<-c("-")
  # permet de synthétiser la p.value sous forme d'un code à étoiles

  value=c(export[i,4],export[i,2],export[i,3],export[i,1])
  group=c(paste("Maintien: ",export[i,4]," mailles",sep=''),
          paste("Apparition: ",export[i,2]," mailles",sep=''),
          paste("Disparition: ",export[i,3]," mailles",sep=''),
          paste("Absence: ",export[i,1]," mailles",sep=''))
  Tab=data.frame(group,value)
  Titre<-NULL
  if (Applicable[i] == 0) Titre<- c("Non Applicable") else
    Titre<- paste(round(mcnemar.test(mat)$statistic,3)," (",sign,")",sep='')
  photo<-paste(getwd(),rownames(export[i,]),".jpg",sep='')

colors <- c("#FFFF00","#FF0000","#0000FF","#FF00FF")
group=c("Maintien","Apparition","Disparition","Absence")

  jpeg(filename=photo,res=600,width=6,height=6,units="cm")
  treemap(Tab,
          index="group",
          vSize="value",
          type="index",
          palette=colors,
          fontsize.labels = 5,
          force.print.labels = T,
          fontface.labels = 4,
          title= Titre,
          fontsize.title = 10,
          border.lwds= 0.5,
          fontcolor.labels="black")
  dev.off()

  Linkimg[i]<-photo
}

export1<-cbind(tabl,Somme,stat_mcnemar,pvalue_mcnemar,Linkimg)
write.csv(export1,paste(getwd(),"/res-mcnemar.csv",sep=''))
0/5 - [0 rating]




Cartographie R [R]

27.08.2018     Remy_Moine      cartographie points polygones vecteurs SIG 

  Un essai de synthèse des éléments rassemblés pour obtenir une carte de différentes données vecteur sous R.

Il est basé en partie sur les éléments donnés par la R-Graph Gallery (https://www.r-graph-gallery.com/). De plus amples explications sont disponibles sur ce site.

Les données en entrée doivent suivre une projection Mercator (EPSG 3857).

packages additionels: "ggplot2,ggmap, broom, sp, rgdal, maptools, plyr, rgeos"

R version 3.4.4
OS: Linux Xubuntu 16.04
library(ggmap)
library(ggplot2)
library(broom)
library(sp)
library(rgdal)
library(maptools)
library(plyr)
library(rgeos)

# stockage des projections usuelles
l93a<-"+init=epsg:27572"
l93<-"+init=epsg:2154"
wgs84<-"+init=epsg:4326"
merca<-"+init=epsg:3857"

# Chargement d'une couche de polygones représentant des secteurs de prospection
secteurs<-readOGR( dsn= getwd() , layer="secteurs") 
proj4string(secteurs)<-CRS(l93)
# Extraction des vertex du polygone pour être représenté sous ggplot.
secteurs@data$id = rownames(secteurs@data)
secteurs.points = fortify(secteurs,region="id")
secteurs.df = join(secteurs.points, secteurs@data,by="id")

# Création d'un point géolocalisé d'inventaire
# !ATTENTION ! Tout doit être en epsg 3857
centro<-as.data.frame(cbind(lon=c(44.95),lat=c(6.10)))# tableau contenant les points EN MERCATOR (EPSG=3857 !!!)

# Récupération d'un fond cartographique googlemap satellite
# Le fond cartographique couvre la zone spécifiée entre guillemet (donc ici la France). La recherche s'effectue comme dans Googl...aps.
map <- ggmap(get_googlemap("France", zoom = 10, maptype = "satellite"))# zoom à 0-> terre entière, zoom à 20-> champs en détail, 

# Assemblage final
carte<-map +
  geom_point(data=centro,aes(x=lat,y=lon,colour=rownames(centro),shape=rownames(centro)),size=3)+
  xlab("Longitude")+ ylab("Latitude")+
  scale_color_manual(values=c("blue"))+
  theme(legend.position='none')
0/5 - [0 rating]




Lister les screens ouverts sur des serveurs distants [Bash]

28.08.2018     Kumquatum      tool bashrc screen 

  Affiche le nom de vos screens ouvert sur une liste de serveurs donnés
À ajouter dans le .bashrc pour une utilisation rapide
function screenAllServ ()
{
echo -e "#######################\n### Current screens ###\n#######################"
for servName in Serv1 Serv2 Serv3
do
ssh username@$servName.adresse.de.mon.serveur.com "echo -e \"----\n<><><> \"\$(hostname)\"<><><>\"; ls /var/run/screen/S-username/ | cat"
done
} 
0/5 - [0 rating]




Simple Fast Kmer counter (with small k) [C++]

28.01.2019     Natir      cpp kmer counter 

  A simple and fast kmer counter.

limitation:
- k must always be odd
- k must be choose at compile time
- the memory usage is 2 ^ (k*2-1) bytes
- in this version the maximum count of a kmer is 255 (see comments)
- in this version the maximum value of k is 31 (see comments)
- forward and reverse kmers are counted at the same time


The main function provides a simple example.
#include <bitset>

template <size_t K>
using kmer_t = typename std::conditional<
  K <= 4 && K % 2 == 1,
  std::uint_least8_t,
  typename std::conditional<
    K <= 8 && K % 2 == 1,
    std::uint16_t,
    typename std::conditional<
      K <= 16 && K % 2 == 1,
      std::uint32_t,
      typename std::conditional<
        K <= 32 && K % 2 == 1,
        std::uint64_t,
        std::false_type
      >::type
    >::type
  >::type
>::type;

template <std::uint8_t K>
constexpr std::uint64_t comp_mask() {
  if(K <= 4 && K % 2 == 1) {
    return 0b10101010;
  } else if(K <= 8 && K % 2 == 1) {
    return 0b1010101010101010;
  } else if(K <= 16 && K % 2 == 1) {
    return 0b10101010101010101010101010101010;
  } else if(K <= 32 && K % 2 == 1) {
    return 0b1010101010101010101010101010101010101010101010101010101010101010;
  } else {
    return 0;
  }
}

template <std::uint8_t K>
constexpr std::uint64_t max_k() {
  if(K <= 4) {
    return 4;
  } else if(K <= 8) {
    return 8;
  } else if(K <= 16) {
    return 16;
  } else if(K <= 32) {
    return 32;
  } else {
    return 0;
  }
}	     

template<std::uint8_t K>
kmer_t<K> seq2bit(std::string seq) {
  kmer_t<K> kmer = 0;
  
  for(auto n: seq) {
    kmer <<= 2;
    kmer |= ((n >> 1) & 0b11);
  }
  
  return kmer;
}

template<std::uint8_t K>
kmer_t<K> reverse_2(kmer_t<K> kmer) {
  kmer_t<K> reversed = 0;

  for(kmer_t<K> i = 0; i < K-1; ++i) {
    reversed = (reversed ^ (kmer & 0b11)) << 2;
    kmer >>= 2;
  }
  
  return reversed ^ (kmer & 0b11);
}

template<std::uint8_t K>
uint8_t parity(kmer_t<K> kmer) {
  return std::bitset<K*2>(kmer).count() % 2; 
}
 
template<std::uint8_t K>
kmer_t<K> get_cannonical(kmer_t<K> kmer) {
  uint8_t cleaning_move = (max_k<K>() - K) * 2;

  if(parity<K>(kmer) == 0) {
    return kmer;
  } else {
    
    kmer_t<K> hash = (reverse_2<K>(kmer) ^ comp_mask<K>());
    hash <<= cleaning_move;
    hash >>= cleaning_move;
        
    return hash;
  }
}

template<std::uint8_t K>
kmer_t<K> hash_kmer(kmer_t<K> kmer) {
  return get_cannonical<K>(kmer) >> 1;
}

#include <iostream>
#include <vector>

int main(int argc, char* argv[]) {

  constexpr std::uint8_t k = 3; // must be odd and minus than 32 change test at line 119

  std::string seq = "AAACCCTTTGGG";

  // To increase maximal count change uint8_t, to unint16_t, uint32_t or uint64_t
  std::vector<std::uint8_t> kmer2count(1 << (k * 2 -1), 0); // 1 << n <-> 2^n

  for(size_t i = 0; i != seq.length() - k + 1; i++) {
    std::string sub_seq = seq.substr(i, k);
    
    uint32_t hash = hash_kmer<k>(seq2bit<k>(sub_seq));

    if(kmer2count[hash] != 255) {
      kmer2count[hash]++;
    }
  }
  
  for(size_t i = 0; i != kmer2count.size(); i++) {
    std::cout<<int(kmer2count[i])<<" ";
  }
  std::cout<<std::endl;
  
  return 0;
}
5/5 - [1 rating]




Average timeit [Python]

26.02.2019     Yo_O      timer python chronophage temps benchmark 

  Une fonction timeit permettant de placer @timeit en decorateur d'une methode à chronométrer.
Petit bonus : joue n fois la méthode décorée et rend le temps moyen.

source originale : https://github.com/realpython/materials/blob/master/pandas-fast-flexible-intuitive/tutorial/timer.py
import functools
import gc
import itertools
import sys
from timeit import default_timer as _timer


def timeit(_func=None, *, repeat=3, number=1000, file=sys.stdout):
    """Decorator: prints time from best of `repeat` trials.
    Mimics `timeit.repeat()`, but avg. time is printed.
    Returns function result and prints time.
    You can decorate with or without parentheses, as in
    Python's @dataclass class decorator.
    kwargs are passed to `print()`.
    >>> @timeit
    ... def f():
    ...     return "-".join(str(n) for n in range(100))
    ...
    >>> @timeit(number=100000)
    ... def g():
    ...     return "-".join(str(n) for n in range(10))
    ...
    """

    _repeat = functools.partial(itertools.repeat, None)

    def wrap(func):
        @functools.wraps(func)
        def _timeit(*args, **kwargs):
            # Temporarily turn off garbage collection during the timing.
            # Makes independent timings more comparable.
            # If it was originally enabled, switch it back on afterwards.
            gcold = gc.isenabled()
            gc.disable()

            try:
                # Outer loop - the number of repeats.
                trials = []
                for _ in _repeat(repeat):
                    # Inner loop - the number of calls within each repeat.
                    total = 0
                    for _ in _repeat(number):
                        start = _timer()
                        result = func(*args, **kwargs)
                        end = _timer()
                        total += end - start
                    trials.append(total)

                # We want the *average time* from the *best* trial.
                # For more on this methodology, see the docs for
                # Python's `timeit` module.
                #
                # "In a typical case, the lowest value gives a lower bound
                # for how fast your machine can run the given code snippet;
                # higher values in the result vector are typically not
                # caused by variability in Python’s speed, but by other
                # processes interfering with your timing accuracy."
                best = min(trials) / number
                print(
                    "Best of {} trials with {} function"
                    " calls per trial:".format(repeat, number)
                )
                print(
                    "Function `{}` ran in average"
                    " of {:0.3f} seconds.".format(func.__name__, best),
                    end="\n\n",
                    file=file,
                )
            finally:
                if gcold:
                    gc.enable()
            # Result is returned *only once*
            return result

        return _timeit

    # Syntax trick from Python @dataclass
    if _func is None:
        return wrap
    else:
        return wrap(_func)
5/5 - [1 rating]




recherche d une sequence mystere a l aide de biopython [Python]

04.04.2019     erwan06      biopython fasta arabette photosynthese 

  ce snippet reprend le code de recherche de la protéine codée vraisemblable comme cela est décrit sur le site http://arn16s.ovh (étape n°3) à partir de la séquence fasta "U91966". La protéine codée par l'une des six phases est alors affichée (Rubisco)
import Bio
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.Alphabet import generic_dna

#on importe la sequence mystere avec les outils en ligne

def seq_codante_de_phase(phase):
    meilleure_seq=compteur=i=j=start=stop=0
    longueur = len(phase)-1
    for compteur in range (i,longueur):
        if (phase[i] != 'M'):
            i += 1
        else:
            start = i
            j = start
            while (phase[j] != '*' and j < len(phase)-1):
                j += 1
            if (phase[j] == "*"):
                stop = j
            seq_codante = phase[i:j]
            if (len(seq_codante) > meilleure_seq):
                meilleure_seq = len(seq_codante)
                meilleur_start = start
                meilleur_stop = stop
                proteine_de_phase = phase[meilleur_start+1:meilleur_stop]
            i += 1
    return(proteine_de_phase)        

#main

mysterious_sequence = SeqIO.read(open('my_sequence.fasta'),
'fasta',
alphabet=generic_dna).seq


phase_1 = mysterious_sequence[0::]
phase_2 = mysterious_sequence[1::]
phase_3 = mysterious_sequence[2::]

#on renverse la sequence mystere pour creer les 3 dernieres phases 
complement_sequence = mysterious_sequence.complement()
reverse_sequence = complement_sequence[::-1]

#sequence à partir du dernier nucleotide en sens contraire
phase_4 = reverse_sequence[0::]
#sequence à partir de l'avant-dernier nucleotide en sens contraire
phase_5 = reverse_sequence[1::]
#sequence à partir de l'antépénultième nucleotide en sens contraire 
phase_6 = reverse_sequence[2::]

prot_1 = str(phase_1.translate())
prot_2 = str(phase_2.translate())
prot_3 = str(phase_3.translate())
prot_4 = str(phase_4.translate())
prot_5 = str(phase_5.translate())
prot_6 = str(phase_6.translate())

liste=[seq_codante_de_phase(prot_1),seq_codante_de_phase(prot_2),seq_codante_de_phase(prot_3),seq_codante_de_phase(prot_3),seq_codante_de_phase(prot_4),seq_codante_de_phase(prot_5),seq_codante_de_phase(prot_6)]
sorted(liste, key=len)
#la phase vraisemblable est la plus longue, donc la première de la liste triee par longueur
print (liste[1])
0/5 - [0 rating]




Couleur moyenne d'un vecteur de couleur en RGB + alpha [R]

30.07.2019     gdevailly      R RGB HTML couleurs 

  Petite fonction retournant la couleur moyenne d'un vecteur de couleur en espace RGB (rouge, vert, bleu) + canal alpha (transparence).
Il y a des espaces colorimétriques sans plus adapté à ce genre d'opérations.
mean_color <- function(mycolors) {
    R     <- strtoi(x = substr(mycolors,2,3), base = 16)
    G     <- strtoi(x = substr(mycolors,4,5), base = 16)
    B     <- strtoi(x = substr(mycolors,6,7), base = 16)
    alpha <- strtoi(x = substr(mycolors,8,9), base = 16)

    return(
        rgb(
            red   = round(mean(R)),
            green = round(mean(G)),
            blue  = round(mean(B)),
            alpha = round(mean(alpha)),
            maxColorValue = 255
        )
    )
}

mean_color(c("#000000FF", "#FFFFFFFF"))
mean_color(c("#FF0000FF", "#00FF00FF", "#FFFF00FF", "#FF0000FF"))
mean_color(rainbow(8))
0/5 - [0 rating]




exemple d'utilisation du module "pyensembl" via un NoteBook Jupyter (Anaconda) [Python]

25.09.2019     erwan06      Anaconda Jupyter pyensembl 

  on va comparer les références de transcrits (identifiants Ensembl grch38) avec le fichier de référence "VariantNGS.tsv" suivant:
Gene name Accession Number
Gene name Accession Number
VHL ENST00000256474
SETD2 ENST00000409792
PBRM1 ENST00000337303
FHIT ENST00000341848
RASSF1 ENST00000357043
KDM5C ENST00000375401
MITF ENST00000394351
BAP1 ENST00000460680
KDM6A ENST00000377967
ABL1 ENST00000318560
AKT1 ENST00000349310
ALK ENST00000389048
APC ENST00000457016
ATM ENST00000278616
BRAF ENST00000288602
CDH1 ENST00000261769
CDKN2A ENST00000304494
CSF1R ENST00000286301
CTNNB1 ENST00000349496
DDR2 ENST00000367922
EGFR ENST00000275493
ERBB2 ENST00000269571
ERBB4 ENST00000342788
EZH2 ENST00000320356
FBXW7 ENST00000281708
FGFR1 ENST00000447712
FGFR2 ENST00000358487
FGFR3 ENST00000440486
FLT3 ENST00000241453
GNA11 ENST00000078429
GNAQ ENST00000286548
GNAS ENST00000371085
HNF1A ENST00000257555
HRAS ENST00000397596
IDH1 ENST00000345146
IDH2 ENST00000330062
JAK2 ENST00000381652
JAK3 ENST00000458235
KDR ENST00000263923
KIT ENST00000288135
KRAS ENST00000311936
MAP2K1 ENST00000307102
MET ENST00000318493
MLH1 ENST00000231790
MPL ENST00000372470
NOTCH1 ENST00000277541
NPM1 ENST00000517671
NRAS ENST00000369535
PDGFRA ENST00000257290
PIK3CA ENST00000263967
PTEN ENST00000371953
PTPN11 ENST00000351677
RB1 ENST00000267163
RET ENST00000355710
SMAD4 ENST00000342988
SMARCB1 ENST00000263121
SMO ENST00000249373
SRC ENST00000445403
STK11 ENST00000326873
TP53 ENST00000269305


#!/usr/bin/env python
# coding: utf-8

# # Installing PyEnsembl
# `PyEnsembl`is a Python interface to Ensembl reference genome metadata such as exons and transcripts. `PyEnsembl` downloads GTF and FASTA files from the Ensembl FTP server and loads them into a local database. PyEnsembl can also work with custom reference data specified using user-supplied GTF and FASTA files.

# In[52]:

import pyensembl


# In[53]:

from pyensembl import EnsemblRelease


# In[54]:

data = EnsemblRelease(78)


# In[55]:

import os, sys


# In[56]:

os.path.isfile("TruSeq.txt")


# In[57]:

f = open("TruSeq.txt",'r',encoding='utf8')


# In[58]:

L = f.readlines()


# In[59]:

f.close()


# In[60]:

genes = [str(L[i]).rstrip("\n") for i in range(1,len(L))]


# In[61]:

liste = [str(genes[i]).split('\t') for i in range(len(genes))]


# In[62]:

liste_simple = []
for elt in range(0,len(liste)):
    liste_simple = liste_simple + liste[elt]


# In[63]:

print(liste_simple)


# ## get all Genes ID from TrueSeq Amplicon cancer panel

# In[64]:

f = open("TruSeq_geneIDs.txt",'w',encoding = 'utf8')
tTruSeq = []

for i in range(0,len(liste_simple)):
    gene_courant = str(liste_simple[i]).rstrip()
    gene_ID = pyensembl.ensembl_grch38.genes_by_name(gene_courant)[0].gene_id
    f.write("gene: {} , gene_ID : {} \n".format(str(liste_simple[i]),gene_ID))
    tTruSeq = [gene_ID, pyensembl.ensembl_grch38.genes_by_name(gene_courant)[0].transcripts] + tTruSeq
f.close()
    

# In[65]:

f = open("genotype_transcripts.txt",'w',encoding = 'utf8')
f.write(str(tTruSeq))
f.close()


# ## compare genotype Transcript IDs vs NGS phenotypic transcripts (VariantsNGS.tsv)

# In[66]:

import os, sys


# In[67]:

import pyensembl


# In[68]:

os.path.isfile("VariantsNGS.tsv")


# In[69]:

f = open("VariantsNGS.tsv",'r',encoding='utf8')


# In[70]:

L2 = f.readlines() 


# In[71]:

f.close()


# In[72]:

genes = [str(L2[i]).rstrip("\n") for i in range(1,len(L2))]


# In[73]:

liste2 = [str(genes[i]).split('\t') for i in range(len(genes))]


# In[74]:

liste_simple2 = []
for elt in range(0,len(liste2)):
    liste_simple2 = liste_simple2 + liste2[elt]

# In[75]:

print(liste_simple2)

# ## manual cross-check of each transcript IDs

# In[76]:

VHL = pyensembl.ensembl_grch38.genes_by_name("VHL")


# In[77]:

tVHL = VHL[0].transcripts


# In[78]:

print(tVHL[0].transcript_id)


# same Transcript IDs found

# In[79]:

SETD2 = pyensembl.ensembl_grch38.genes_by_name("SETD2")
tSETD2 = SETD2[0].transcripts
print(tSETD2[0].transcript_id)


# different transcript IDs found : ENST00000431180 != ENST00000409792

# l'automatisation de la comparaison des deux listes simples n'est pas encore terminee




0/5 - [0 rating]