staging.inyokaproject.org

Moin, ich muss wöchentlich bei etwa 50 Files die Menge an A T C G und Ns finden und in eine Tabelle eintragen. Eine Arbeit, die auf die Dauer sehr eintönig und frustierend ist.

Nun würde ich gerne ein sh Skript erstellen, das aus allen Files (sogenannte fasta Files) die Info herauszieht wie viele A T C G und Ns zu finden sind. Als Output wünsche ich mir (meinetwegen in der Shell -muss das eh dann in eine Excel-Tabella manuell übetragen) eine Liste mit Probennummer (also File-Name) und dahinter die Anzahl A T G C und Ns.

Also in etwa so:

AB-2021-30-5: A 3500 C 3687 G 4789 T 5678 N 120 AB-2021-30-6: A 3100 C 3908 G 4590 T 5660 N 124

ich habe mal etwas mit grep rumprobiert aber es klappt nicht so wie erwünscht. Gebe ich einen Ordner mit fasta Files an, wird in der Shell die Gesamtmenge an z.B. Ts in allen Files zusammengezählt und angegeben und nicht für jedes File einzeln:

1
2
3
4
5

FILES="/NAS/Auswertung/Auftrag/2021/ABC/Sequenzen/FASTA/*.fasta"

for file in $FILES; do
    grep -o T $FILES | wc -l
done

Ich habe auch noch ein paar andere Dinge ausporbiert, aber komme irgendwie nicht weiter. Vielleicht hat ja jemand von euch eine Idee.

Eine Möglichkeit wäre mit tr.

1
2
3
4
5

echo A $(tr -c -d A < input.txt | wc -c)
echo C $(tr -c -d C < input.txt | wc -c)
echo G $(tr -c -d G < input.txt | wc -c)
echo T $(tr -c -d T < input.txt | wc -c)
echo N $(tr -c -d N < input.txt | wc -c)

Falls die Dateien sehr groß sind solltest du jedoch eine effizientere Lösung finden, die nur einmal die Datei einliest statt für jeden Buchstaben von vorne.

PS: Bei deinem for loop ist komisch, daß du im loop dann mit FILES arbeitest statt mit der Loopvariable file

Hm, warum nutzt Du eine for-Schleife, die Du eigentlich doch nicht nutzt?

So wäre es dann wohl richtig:

FILES="/NAS/Auswertung/Auftrag/2021/ABC/Sequenzen/FASTA/*.fasta"

for file in $FILES; do
    grep -o T $file | wc -l
done

Dafür gibt es BioPerl und BioPython: https://biopython.org/wiki/Seq

Hm, warum nutzt Du eine for-Schleife, die Du eigentlich doch nicht nutzt? So wäre es dann wohl richtig:

Ah mein Fehler! Danke! Jetzt kriege ich so immerhin schonmal die Anzahl an Ts untereiander für jedes File. Leider ohne Namen der Files davor. -H klappt da auch nicht - es wird kein Filename angezeigt im Output

Dafür gibt es BioPerl und BioPython: https://biopython.org/wiki/Seq

Danke! Gucke ich mir mal an

Falls die Dateien sehr groß sind solltest du jedoch eine effizientere Lösung finden, die nur einmal die Datei einliest statt für jeden Buchstaben von vorne.

Nein, groß sind sie eigentlich nicht. Ich probiere das mal aus. Danke!

Hm, wenn das -H bei grep nicht funktioniert, gib den Dateinamen einfach vorher nochmal aus:

FILES="/NAS/Auswertung/Auftrag/2021/ABC/Sequenzen/FASTA/*.fasta"

for file in $FILES; do
    echo -n "$file "; 
    grep -o T $file | wc -l;
done

Ich habe es genauso hinbekommen wie ich es wollte. Bekomme nun eine Liste für jede Probe:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

read -p "ID Auftrag eingeben:" ID 

cd /NAS/Auswertung/Auftrag/2021/ABC/"$ID"/Analysen/Pipe/results/consensuses_iupac

LIST="$(find *.fasta)"

for file in $LIST; do
	echo "$file"
echo A $(tr -c -d A < $file | wc -c)
echo C $(tr -c -d C < $file| wc -c)
echo G $(tr -c -d G < $file| wc -c)
echo T $(tr -c -d T < $file| wc -c)
echo N $(tr -c -d N < $file | wc -c)
echo ""
done

cd

Vom C muss ein C abgezogen werden, weil in der ersten Zeile im Filename ein C vorkommt. Es wäre "nice to have" aber nicht unbedingt nötig, wenn direkt ein C abgezogen werden würde. Da wir es aber so oder so in eine Excel Tabelle auf unseren Windows Rechner übertragen müssen und wir es so gerade eben noch hinbekommen -1 zu rechnen, wäre es nicht nötig.

Auch eine zweite Sache wäre schön, aber nicht unbeding notwendig. Das Skript ist an der Stelle oben noch nicht zu Ende, denn ich brauche auch die Ns eines anderen Files. Diese Anzahl der Ns oben muss von der Anzahl Ns unten abgezogen werden. Ist nicht weiter schwierig im Kopf zu rechnen, denn meistens ist das Ergebnis 1 oder 2. Und wenn es mal ein paar mehr sind, ist es auch kein Problem es auszurechnen. Aber es wäre natürlich schön, wenn man direkt die Zahl da hätte, die man eintragen muss...

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

cd NAS/Auswertung/Auftrag/2021/ABC/"$ID"/Analysen/Pipe/results/consensuses_masked

LIST="$(find *.fasta)"


for file in $LIST; do
	echo "$file"
echo N $(tr -c -d N < $file | wc -c)
echo ""
done
cd

Der Output sieht so aus:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

----------A T C G N consensuses_iupac:----------

Probe_1_.iupac_consensus.fasta
A 7864
C 4824
G 5120
T 8465
N 3622

Probe_2_.iupac_consensus.fasta
A 7047
C 4336
G 4623
T 7670
N 6211

----------N consensuses_masked:----------

Probe_1_.masked_consensus.fasta
N 3630

Probe_2_.masked_consensus.fasta
N 6214

Im Kopf müssten wir also 3630-3622 rechnen. Ist kein Problem 😉 Aber unten wäre es eben prima, wenn da direkt statt 3630 einfach 8 stehen würde.

Bis auf die Reihenfolge (hier: Anzahl vorne, Suchmuster hinten) lässt sich die Aufgabe mit diesem Einzeiler lösen:

grep -H -e '[ATGCN]' -o -r VERZEICHNIS | sort | uniq -c 

Für VERZEICHNIS ist natürlich der Ordner anzugeben, in dem sich die zu durchsuchenden Textdateien gewöhnlich aufhalten.

Vielen Dank kB!!!

es klappt aber Output gefällt mir nicht:

1
2
3
4
5

8775 Probe_.iupac_consensus.fasta:A
5395 Probe_.iupac_consensus.fasta:C
5775 Probe_.iupac_consensus.fasta:G
510 Probe_.iupac_consensus.fasta:N
9435 Probe_.iupac_consensus.fasta:T

da ist das unten wesentlich übersichtlicher

Wegen Rechnen schau Dir mal [Shell/Bash-Skripting-Guide_für_Anfänger/#Rechnen-mit-der-Shell:] an.

Wenn das ATCGN nur ein bestimmter Bereich in der Datei ist (z.B. erst ab Zeile zwei) dann kannst das z.B. mit tail -n +2 direkt so ohne die störende Zeile ausgeben. Statt da hinterher was nachzurechnen.

Das tr | wc ist ein Stream, das grep | sort macht erstmal alles im Arbeitsspeicher. Wenns nur um ein paar tausend Zeilen geht, wirds wohl nicht ins Gewicht fallen. Wenn du so ne Datei hast mit paar Gigabyte dann funktioniert der grep sort Ansatz irgendwann nicht mehr.

Am Ende sind das paar Zeilen C / Python um das ordentlich in einem Rutsch zu zählen. Irgendwann sollte man das einfach nicht mehr in der Shell machen.

SK-N-BE schrieb:

Vielen Dank kB!!!

es klappt aber Output gefällt mir nicht:

Die Formatierung kannst du dir doch selbst basteln, wie du willst. Also z.B.

grep -H -e '[ATGCN]' -o -r VERZEICHNIS | sort | uniq -c | gawk -F"[ :]" 'BEGIN{f="File";split("A C G N T",a)
for(i in a){arr[a[i]]=a[i]}}
function p(){printf "%s;",f;for(i in a){printf "%4s;",arr[a[i]]};printf "\n"}
{if($2==f){arr[$3]=$1}else{p();delete arr;f=$2;arr[$3]=$1}
}END{p()}' | column -s";" -t
#
# Beispiel:
#
echo "8775 Probe_.iupac_consensus.fasta:A
5395 Probe_.iupac_consensus.fasta:C
5775 Probe_.iupac_consensus.fasta:G
510 Probe_.iupac_consensus.fasta:N
9435 Probe_.iupac_consensus.fasta:T
7864 Probe_1_.iupac_consensus.fasta:A
4824 Probe_1_.iupac_consensus.fasta:C
5120 Probe_1_.iupac_consensus.fasta:G
8465 Probe_1_.iupac_consensus.fasta:T
3622 Probe_1_.iupac_consensus.fasta:N
7047 Probe_2_.iupac_consensus.fasta:A
4336 Probe_2_.iupac_consensus.fasta:C
4623 Probe_2_.iupac_consensus.fasta:G
7670 Probe_2_.iupac_consensus.fasta:T
6211 Probe_2_.iupac_consensus.fasta:N" |
gawk -F"[ :]" 'BEGIN{f="File";split("A C G N T",a)
for(i in a){arr[a[i]]=a[i]}}
function p(){printf "%s;",f;for(i in a){printf "%4s;",arr[a[i]]};printf "\n"}
{if($2==f){arr[$3]=$1}else{p();delete arr;f=$2;arr[$3]=$1}
}END{p()}' | column -s";" -t
File                               A     C     G     N     T
Probe_.iupac_consensus.fasta    8775  5395  5775   510  9435
Probe_1_.iupac_consensus.fasta  7864  4824  5120  3622  8465
Probe_2_.iupac_consensus.fasta  7047  4336  4623  6211  7670

Ein alternativer Ansatz:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

#!/usr/bin/ruby

ARGV.each do |file|
  counts = Hash.new 0

  File.foreach file do |line|
    line.scan /[ACGTN]/ do |ch|
      counts[ch] += 1
    end
  end

  printf '%s:', file

  %w{A C G T N}.each do |ch|
    printf ' %s %d', ch, counts[ch]
  end

  puts
end

Zu Erklären bin ich gerade zu müde - sorry. Vorteil so einer Skriptsprache ist, dass man damit jede Menge Dinge sehr flexibel lösen kann. Ruby hat vor allem den Vorteil, dass reguläre Ausdrücke gut in die Sprache integriert sind. Wenn Du viel in Dateien suchen musst, wirst Du vermutlich früher oder später um Regex nicht herum kommen. 😉

Vielen Dank für eure tolle Hilfe!!!

Shiro:

vielen Dank! Ich verstehe den Code nicht. Habe ihn aber trotzdem probiert und der Rechner rechnet und rechnet und rechnet ohne dass ein Output angegeben wird. Mit dem, was ich bisher gemacht habe, bekomme ich die Ergebnisse sofort. Schade, dass es nicht klappt, denn der Output - also dass alle Basen nebeneinader stehen- würde mir besser gefallen - dann kann ich das alles nämlich direkt in eine txt Datei überführen, die ich dann mit Exel öffne und die Zahlen direkt in die Excel-Datei kopieren.

Leider habe ich keine Zeit mich noch näher damit beschäftigen. Wenn ich mal Zeit zwischendurch habe, werde ich da aber weiter machen.

rklm: dir auch vielen vielen Dank!!! leider habe ich, wie gerade geschrieben keine Zeit mich näher mit allem zu beschäftigen. Ist aber auf jeden Fallm ein total spannender Ansatz und da ich noch viele viele Jahre hier hoffentlich arbeiten werde und viele spannende Sachen programmieren darf, hoffe ich, dass ich das irgendwann lerne.

Also, ich habe mich nun für die Lösung entschieden, die für unsere Zwecke nun gut funktioniert:

1
2
3
4
5

echo A $(tr -c -d A < $file | wc -c) 
echo C $(tr -c -d C < $file | wc -c)
echo G $(tr -c -d G < $file| wc -c)
echo T $(tr -c -d T < $file| wc -c)
echo N $(tr -c -d N < $file | wc -c)

Wenn das ATCGN nur ein bestimmter Bereich in der Datei ist (z.B. erst ab Zeile zwei) dann kannst das z.B. mit tail -n +2 direkt so ohne die störende Zeile ausgeben. Statt da hinterher was nachzurechnen.

Frostschutz: Vielen Dank für den Tipp. Leider habe ich das nicht hinbekommen. Wo müsste ich denn tail -n +2 einfügen? Wenn ich es so mache, wird in der Zeile 0 angezeigt:

1

echo C $(tr -c -d C | -n +2  < $file | wc -c)

Wenn ich es nach ganz hinten verschiebe, wird für C garnix angezeigt.

Zum Schluss noch einmal Vielen Dank an alle!!! Ihr seid super!