{ Les snippets de bioinfo-fr.net }


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 vote]




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 vote]




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 vote]




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 vote]




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 vote]




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 vote]




Retire le numéro de version du code Ensembl des gènes / transcrits [R]

05.07.2018     gdevailly      GENCODE Ensembl 

  Les codes Ensembl des gènes / transcrits viennent parfois avec un chiffre final traçant le nombre de modifications de ces annotations:

Par exemple, l'annotation de ENSG00000134046.11 a été modifié 11 fois.

Cette petit expression régulière permet de s'en débarrasser, pour faciliter les croisements de listes de gènes.
x <- c("ENSG00000134046.11", "ENSG00000229807.10")

sub("\\.[0-9]*$", "", x)
# [1] "ENSG00000134046" "ENSG00000229807"
0/5 - [0 vote]




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 votes]




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 votes]




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 vote]





<<< Snippets récents -- Ancients snippets >>>
Powered by canSnippet CE - by ademcan