Script pour enlever des éléments transposable chevauchant et avec un écart de 1000 bp

[Oceane666]

Bonjour à tous,
J'aurais besoin d'aide pour écrire un script qui enlève tous les éléments transposable qui overlaps (qui forme des clusters), et ne garder que les éléments transposable qui n'overlap pas avec d'autre éléments transposable et qui ont une distance minimum en amont et aval d'un autre élément transposable de 1000 bp .

J'ai un fichier avec ce nom : filtered_hypo_Gd293_Dorado_mod basecalling.fasta.out
Et ce type de données :
1739 3.3 0.0 0.0 contig_1 21700 21913 (1813093) + RM2_rnd-1_family-0 LTR/Gypsy 1666 1879 (4848) 1
9413 6.9 4.7 1.2 contig_1 21994 23389 (1811617) C RM2_rnd-1_family-0 LTR/Gypsy (5064) 1663 220 2
1531 1.6 0.0 0.0 contig_1 23260 23444 (1811562) + EDTA_TE_00000254_inc DNA 206 390 (6) 3 *

Les colonnes qui m'intéressent sont la 5,6, et 7
La colonne 5 correspond au nom de mon contig, la 6 à la position de départ de mon élément transposable et la 7 à la position de fin.
Ainsi pour mes 3 lignes, on voit que sur la deuxième ligne mon élément transposable n'overlap pas avec celui de la première ligne mais que l'écart est inférieur à 1000 bp (entre la position end de ma ligne 1 et start de ma ligne 2) je dois donc l'enlever, mais avant de l'enlever je dois aussi comparer cette ligne à la ligne du dessous (donc la ligne 2 avec la 3) , ici je vois un chevauchement entre la ligne 3 et 2 je dois donc aussi l'enlever. Ainsi dans mes 3 lignes je n'en garderai ici aucune car elle n'ont pas un écart minimum de 1000 bp en amont et aval et possède des overlap.

Merci d'avance pour votre aide !
-----

[Henrix]

Bonjour,
Sur ce fil-->: https://forums.futura-sciences.com/p...51-lignes.html
On t'a donné trois solutions.
Tu n'a pas daigner faire retour (répondre).
Si tu es sous Windows, je peux te faire un soft dédié, open source et sans installation, qui d'ailleurs pourra intégrer tes trois demandes ( https://forums.futura-sciences.com/p...algorithmique/ ).
Mais si tu ne réponds pas... Tu va démotiver les plus courageux.

[pm42]

Tu n'a pas daigner faire retour (répondre).

C'est faux : il a répondu très poliment.

Si tu es sous Windows, je peux te faire un soft dédié

Il utilise bash donc il est peut-être sous Windows mais la probabilité est plus basse.
Et pourquoi faire un soft dédié quand ce qu'il demande s'écrit en général en quelques lignes de script ?

[pm42]

La colonne 5 correspond au nom de mon contig, la 6 à la position de départ de mon élément transposable et la 7 à la position de fin.
Ainsi pour mes 3 lignes, on voit que sur la deuxième ligne mon élément transposable n'overlap pas avec celui de la première ligne mais que l'écart est inférieur à 1000 bp (entre la position end de ma ligne 1 et start de ma ligne 2) je dois donc l'enlever, mais avant de l'enlever je dois aussi comparer cette ligne à la ligne du dessous (donc la ligne 2 avec la 3) , ici je vois un chevauchement entre la ligne 3 et 2 je dois donc aussi l'enlever. Ainsi dans mes 3 lignes je n'en garderai ici aucune car elle n'ont pas un écart minimum de 1000 bp en amont et aval et possède des overlap.

C'est relativement facile à faire mais il y a 2 questions :

- que fais tu pour la 1ère et la dernière ligne ?
- est ce que tu ne penses pas que tu devrais apprendre à programmer au lieu de demander au forum de faire ton boulot à ta place quelque part ?

[A voir en vidéo sur Futura]

[polo974]

En perl, en perl, en perl...

Une partie du sommaire du livre "Minimal Perl For UNIX and Linux People"

3 Perl as a (better) grep command
4 Perl as a (better) sed command
5 Perl as a (better) awk command
6 Perl as a (better) find command

[Henrix]

@: pm42 (#3):
Un soft permet d'avoir directement des résultats et d'ajouter des fonctions selon les besoins.
Et plus besoin de taper des lignes (avec les éventuelles fautes de frappe).

[pm42]

@: pm42 (#3):
Un soft permet d'avoir directement des résultats et d'ajouter des fonctions selon les besoins.

C'est dingue. Et on ne peut pas faire ça dans un script ?

Et plus besoin de taper des lignes (avec les éventuelles fautes de frappe).

Tu n'as jamais entendu parler d'un truc qui s'appelle le copier-coller ?

Au fait, tu devrais t'excuser auprès d'Oceane666 pour l'avoir accusé d'impolitesse à tort.

[pm42]

En perl, en perl, en perl...

A la louche :

#!/usr/bin/perl


$line1=undef;
$line2=undef;
while(my $line3 = <STDIN>) {
     if(defined($line1) && defined($line2))
     { 
	 my @parts1 = split(" ", $line1);
	 my @parts2 = split(" ", $line2);
	 my @parts3 = split(" ", $line3);
	 my $end1 = $parts1[6];
	 my $start2 = $parts2[5];
	 my $end2 = $parts2[6];
	 my $start3 = $parts3[5];
	 if($start2 >= $end1+1000 && $end2 <= $start3 - 1000) {
	     print($line2);
	 }
     }
     elsif(defined($line2))
     { 
	 my @parts2 = split(" ", $line2);
	 my @parts3 = split(" ", $line3);
	 my $start2 = $parts2[5];
	 my $end2 = $parts2[6];
	 my $start3 = $parts3[5];
	 if($end2 <= $start3 - 1000) {
	     print($line2);
	 }
     }
     $line1 = $line2;
     $line2 = $line3;     
}
my @parts1 = split(" ", $line1);
my @parts2 = split(" ", $line2);
my $end1 = $parts1[6];
my $start2 = $parts2[5];
if($end1 <= $start2 - 1000) {
     print($line2);
}

On doit pouvoir faire nettement mieux mais je suis rouillé en perl.

[Oceane666]

Merci infiniment pour votre aide. J'essaye de comprendre le script pm42, mais j'avoue que j'ai un peu de mal sur certaines lignes :
while (my $line3 = <$fh>), sur cette ligne à quoi correspond fh en fait ?
my @parts2 = split(" ", $line2); et quel est la fonction de split exactement, quand l'utiliser

Encore merci pour votre aide, je vais regarder plus en détail des guide de programmation en perl et bash est ce que vous auriez des recommandations personnelles?

[pm42]

Le script en Perl, c'était juste un clin d'œil à polo974.

Je te déconseille d'investir dans ce langage qui est daté. Le mieux de nos jours pour commencer à programmer et manipuler des données c'est dans doute Python.

[pm42]

Et en python, c'est plus concis et rigolo même si on peut là aussi faire mieux :

import pandas as pd
import sys

df = pd.read_csv(sys.stdin, sep=' ', header=None,  names=range(16))

df = df[(df[5] >= df.shift(1)[6]+1000) & (df[6] >= df.shift(-1)[5]-1000)]
df.fillna('', inplace=True)
df = df.apply(lambda x: ' '.join([str(y) for y in x]), axis=1)
for i in range(len(df)):
    print(df.iloc[i])

[Flyingbike]

Le script en Perl, c'était juste un clin d'œil à polo974.

Je te déconseille d'investir dans ce langage qui est daté. Le mieux de nos jours pour commencer à programmer et manipuler des données c'est dans doute Python.

de toute façon c'est ce qui est utilisé pour ce genre de manipulation, où ou trouve une kilochiée d'outils déjà existants et une grosse communauté.

[polo974]

Oui, le coup du perl, c'était de la blague, (jamais je ne coderai en perl) et un peu de provoc, car là il n'y a que des demandes de solutions toutes cuites.

Bref, python est un bon langage pour faire un peu de tout (et le reste, voire n'importe quoi).
Attention, pandas est une lib qui s'appuie sur numpy, où on pense matrice, ce qui permet une écriture très simple sans se palucher le parcours en long et en large des éléments de la (ou des) matrice(s).

Au fait, c'est quoi un bp ???

Sinon, les 2 dernières lignes peuvent être remplacées par:

for lig in df.iloc: print(lig)

[Flyingbike]

bp = base pair = paire de bases

[polo974]

Merci,
une fois dit, ça a du sens...

[pm42]

for lig in df.iloc: print(lig)

En effet.
Ou encore plus court :

df.apply(print)

On peut même faire tout le script en 3 lignes :

import pandas as pd
df = pd.read_csv('foo.dat', sep=' ', header=None,  names=range(16))
df[(df[5] >= df.shift(1)[6]+1000) & (df[6] >= df.shift(-1)[5]-1000)].fillna('').apply(lambda x: print(' '.join([str(y) for y in x])), axis=1)

[polo974]

Trop fort, pandas...
En fait ma remarque était due au fait qu'en general, un truc genre:
for i in range(len(truc)):

Si truc est bien foutu, on peut faire mieux.
J'avoue avoir largué plein de boucles de ce ganre quand j'ai commencé à coder en python (réminiscences du C...).

Sinon, je me demandais, est-il possible d'avoir un itérateur sur le fichier même (avec une fonction "d'oubli" des lignes totalement exploitées et/ou un nombre de lignes à garder en mémoire) tout en ayant l'usage global genre apply?
Le but étant de pouvoir traiter un flux au fil de l'eau ou un giga fichier sans exploser la mémoire.

(some time later...)
Bon, je me suis dis, cherche avant de poster...

Il y a donc le paramètre chunksize sur pandas, mais en plus, il y a une lib "plus mieux" qui s'appuie sur pandas: dask

Il faut vraiment que je creuse ça...

Bon, et maintenant la question qui remet une couche: et pour ouvrir une série de fichiers (genre log journaliers par ex)??? (bon, je saurai créer la fonctionnalité en quelques lignes et par exemple un pipe, mais dans la lib même, ce serait cool aussi)

[pm42]

En fait ma remarque était due au fait qu'en general, un truc genre:
for i in range(len(truc)):
Si truc est bien foutu, on peut faire mieux.

Oui et je n'en fait normalement plus mais dans le cas présent, c'était le plus "intuitif" à coder vite.

Sinon, je me demandais, est-il possible d'avoir un itérateur sur le fichier même (avec une fonction "d'oubli" des lignes totalement exploitées et/ou un nombre de lignes à garder en mémoire) tout en ayant l'usage global genre apply?

Bien sur : tu fais un genérateur avec yield et tu utilises les "for comprehension". Ou la fonction "map".

Le but étant de pouvoir traiter un flux au fil de l'eau ou un giga fichier sans exploser la mémoire.

Sur les OS récents (depuis même quelques décennies), ce n'est pas forcément un sujet : les fichiers ne sont pas lus tels quels mais "mappés" dans la mémoire virtuelle. Donc on peut à la fois le voir en entier sans consommer de la mémoire réelle.

Il y a donc le paramètre chunksize sur pandas, mais en plus, il y a une lib "plus mieux" qui s'appuie sur pandas: dask

dask, c'est surtout pour faire du bigdata.

Bon, et maintenant la question qui remet une couche: et pour ouvrir une série de fichiers (genre log journaliers par ex)??? (bon, je saurai créer la fonctionnalité en quelques lignes et par exemple un pipe, mais dans la lib même, ce serait cool aussi)

Tout dépend de ce que tu veux faire. J'aurais tendance à faire du parallèle soit avec le multiprocess actuel, soit si les I/O prennent le plus gros en en multithreading avec asyncio en attendant les futures versions de Python qui supporteront du vrai multithreading.

Ou alors, j'utiliserais un langage digne de ce nom.