서열자료조작을위한실용적인스크립트예제 작성자 : 정해영 hyjeong@kribb.re.kr jeong0449@gmail.com ( 오류정정이나개선사항에대한의견환영!) 최초작성일 : 2015년 8월 6일최종수정일 : 2015년 9월 22일최초공개일 : 2015년 8월 7일 일러두기 Ÿ 이탤릭체로쓴것은사용자가상황에맞게입력해야하는곳을나타냄 Ÿ 본자료의활용을위해서는리눅스시스템에 EMBOSS 패키지와 samtools가설치되어있어야합니다. Fasta/fastq 파일에몇개의서열이있는지세어보기 $ grep -c '>' multi_sequences.fasta $ awk 'END {print NR/4' file.fastq [awk]: 중괄호앞뒤에공백이있어도되고없어도된다. 죽, {action 이나 { action 어느형태를써도상관이없다. Fastq 파일에서앞부분의 read 4백만개 (1600만줄 ) 만뽑아서별도의파일로저장하기 $ awk 'NR == 1, NR == 16000000' in_file.fastq > out_file.fastq $ head -n 16000000 in_file.fastq > out_file.fastq [awk] 비교연산자앞뒤에공백이없어도된다. Fastq 파일을 (single line) fasta 파일로변환하기 Fasta file 은읽기쉽도록 60 글자마다줄바꿈을하는것이일반적이다. single line fasta 란 fastq 파일처럼줄 바꿈없이서열을하나의줄에표현한것과같다. single line fasta 로전환하면후속작업이편리해진다. $ awk 'NR%4 == 1 {a=substr($0,2); NR%4 == 2 {print ">" a "\n" $0' file.fastq > file_single_line.fasta [awk]: 중괄호 { 내부에명령어가하나만들어가는경우세미콜론을생략해도무방하다. 가독성을높이기위하 여중괄호의바로내부에공백을넣어도된다. [awk]: {print ">" a "\n" $0 과 {print ">"a"\n"$0 은결과가같다. print 명령어로인쇄할인자들 사이에콤마를찍으면 OFS 변수 (= 구분자, 기본값은공백 ) 를사이에넣는다. 만약탭 ( \t ) 을구분자로 설정하고싶으면 awk vofs \t pattern {action 형태로실행하라. $ sed -n '1~4{s/^@/>/;p; 2~4p' file.fastq Multi-line fasta 파일을 single line fasta 파일로변환하기 $ awk '!/^>/ { printf "%s", $0; n = "\n" /^>/ { print n $0; n = "" END { printf "%s", n ' file.fasta > file_single_line.fasta Single line fasta 파일을다시 multi-line fasta 파일로바꾸기 EMBOSS 가설치되어있다면 seqret 유틸리티를이용한다.
$ seqret -sequence single_line.fasta -outseq multi-line.fasta -sequence 와 -outseq 는생략하고입출력파일명만순서에맞게제공해도된다. Fasta 파일 (single line) 내의서열을읽어서 ID 와길이를출력하기 $ awk 'NR%2 == 1 {id=substr($0,2) NR%2 == 0 {print id, " ", length($1)' single_line.fa 길이에따라정렬하여출력하려면 sort 명령어를사용한다. 길이에따른내림차순으로정렬 ( 긴서열이앞으로나 오도록 ) 하려면 sort -k 2 -nr 명령어를실행한다. $ awk 'NR%2 == 1{id=substr($0,2) NR%2 == 0 {print id, " ", length($1)' single_line.fa sort -k 2 -n Fasta 파일 (single line) 읽어서 tab-delimited two-column file(id 와서열 ) 로출력하기 $ awk -vofs="\t" 'NR%2 == 1{id=substr($0,2) NR%2 == 0 {print id, $1' single_line.fa > out.txt ID 와서열로이루어진 two-column file 을 fasta file 로전환하기 $ awk -vofs='' '{print ">",$1,"\n",$2;' out.txt > file.fasta ID 와서열로이루어진 two-column file 에서 1 kb 이상서열만을선택하여 fasta 파일로출력하기 $ awk 'length($2) > 1000 {print ">"$1"\n"$2' file.tab > longerthan1kb.fasta Multi-fasta 파일내의서열들을분리하여각각별도의파일로저장하기 $ awk '/^>/ {OUT=substr($0,2) ".fa"; {print >> OUT; close(out)' multi_sequences.fasta 또는별첨의 fastasplit.pl Perl 스크립트를이용한다. Fasta 파일내에서특정서열만추출하기 (ID 지정 ) 추출할서열의수가많지않아서 one-line script 내에서열 ID 를삽입할수있는경우 $ perl -ne 'if(/^>(\s+)/){$c=grep{/^$1$/qw(id1 id2)print if $c' sequences.fasta 서열 ID 가별도의목록파일에들어있는경우 ( 한줄에 ID 하나 ) $ perl -ne 'if(/^>(\s+)/){$c=$i{$1$c?print:chomp;$i{$_=1 if @ARGV' ids.txt sequences.fasta 1000 bp 이상의서열만추출하여 single line fasta 파일로출력하기 Fasta/fastq 파일을 ID 와서열로구성된 tab-delimited file 로먼저전환해둔다. $ awk 'length($2) > 1000 {print ">"$1"\n"$2' file.tab > longerthan1kb.fasta Fastq 파일로부터 70-120 bp 길이범위의 read 만출력하기
Sequence quality 는고려하지않는단순한방법이므로, quality 를기준으로하는고도한 trimming 이나 filtering 등을위해서는전용도구 ( 예 : FASTX-Toolkit, SolexaQA, prinseq 등 ) 를이용하는것이바람직하다. $ awk 'BEGIN {OFS = "\n" {header = $0 ; getline seq ; getline qheader ; getline qseq ; if (length(seq) >= 70 && length(seq) <= 120) {print header, seq, qheader, qseq' infile.fastq > filtered.fastq 하나의 interleaved fastq file 을 2 개의 paired file 로분리하기 $ awk 'BEGIN{OFS="\n" $0~/\/1$/{header = $0; getline seq; getline qheader; getline qseq; print header, seq, qheader, qseq' infile.fastq > file_1.fastq $ awk 'BEGIN{OFS="\n" $0~/\/2$/{header = $0; getline seq; getline qheader; getline qseq; print header, seq, qheader, qseq' infile.fastq > file_2.fastq 서열의헤더에서 /1 /2와같은 read의 pair 표시가어떻게되어있는지반드시미리확인하여이에맞는 awk 구문내패턴을결정해야한다. @HWI-ST908R:151:D1M9MACXX:1:1101:1373:2469 2:N:0:CCGTCC => $2~/2:/ @BYI6V:00004:00005/1 => $0~/\/1$/ sed 를사용하면한층더간단하게작업할수있다. 서열의 ID 를실제로확인하여패턴을결정해야한다. $ sed -n '/1:N/ {N;N;N;p;' infile.fastq > file_1.fastq $ sed -n '/2:N/ {N;N;N;p;' infile.fastq > file_2.fastq Fasta file 에서특정영역을추출하기 먼저 EMBOSS 의 seqret 을이용하는방법을알아보자. 추출할서열의시작위치는항상끝위치보다앞에있어 야한다. 따라서 reverse strand 에암호화된유전자염기서열을추출할때에도 stop codon 에해당하는위치가 시작위치로주어져야하며, -sreverse 옵션으로상보서열전환을하면된다. $ seqret -sbegin 20 -send 200 single_sequence.fa out.fa $ seqret -sbegin 20 -send 200 single_sequence.fa genbank:out.fa $ seqret -sbegin 20 -send 200 multi_sequences.fa:seq_id out.fa $ seqret -sbegin 20 -send 200 -sreverse single_sequence.fa out.fa seqret 이출력하는서열의 ID 는입력물의것을그대로가져다쓰게된다. 출력서열의 ID 를원하는대로수정하 려면실행시간이좀더소요되지만 sed 를사용하여 > 로시작하는서열 ID 라인을수정하면된다. $ seqret -sbegin 100 -send 200 -stdout -auto multi_sequences.fa:seqid 's/^>.*$/>seqid:100-200/' > out.fa 출력물 : >seqid:100-200 $ seqret -sbegin 100 -send 200 -stdout -auto multi_sequences.fa:seqid 's/^>.\(*\)$/>newid \1/' > out.fa 출력물 : >newid seqid:100-200 < 서열 ID> <start> <end> 형식으로 extract 할서열과위치정보를수록한텍스트파일 (regions) 가존재한다면 저장되어있다면다음과같이하면된다. sed 명령어의치환부분에서변수를큰따옴표와작은따옴표가이중으 로둘러싸고있음에유의하라. $ cat regions while read id start end > do > seqret -stdout -auto -sbegin $start -send $end multi_sequences.fa:$id sed 's/^>.*$/>'"$id:$start-$end"'/' >> out.fa > done 또다른방법으로는 samtools faidx 명령을쓰는것이다. 서열추출전에 fasta 파일에대한인덱스를생성해 두어야한다. EMBOSS 의 seqret 과다른점은출력서열의 ID 는 seqid:start-end 형태로새로씌여진다는것이 다. samtools faidx 로는상보적인서열로전환이되지못하므로필요시에는 EMBOSS 패키지의 revseq 을사용
한다. 필요하다면파일로리다이렉션하기전에위의사례와같이 sed 명령어를사용하여서열의 ID 를수정해도 된다. $ samtools faidx multi_sequence.fa $ samtools faidx multi_sequence.fa seqid 10-200 > out.fa $ revseq out.fa out_rev.fa glimmer 를사용한유전자예측결과물로부터 CDS 서열추출하기 1. NCBI Glimmer site(http://www.ncbi.nlm.nih.gov/genomes/microbes/glimmer_3.cgi) 에세균유전체 서열을업로드하여 glimmer 를실행한뒤결과를텍스트파일로저장하고 fasta 파일과유사하게정리한다. 입력하는서열의 topology 를 circular 로지정하면서열의말단을통과하는유전자구조도예측을하게되므 로주의를요한다. 2. Glimmer 는각 (contig) 서열에대하여동일한 orf 번호를반복하여사용하므로입력물이 multi-fasta 인경우 혼동을유발할수있다. 따라서다음과같은 awk 명령문을이용하여각 orf 에대한유일한식별자를부여하 도록한다. 또한 glimmer 가출력하는유전자의시작부위는번역개시위치기준이다. $ awk 'BEGIN {i=0 /^>/ {id=substr($0, 2) /^orf/ {printf("%s %s ORF%05d %10s %d %d\n", id, $1, i, $2, $3, $4); i++' glimmer_3.txt > glimmer_3.out 3. 유전자위치정보를암호화된 strand 에따라서둘로나누어별도의파일에저장한다. $ awk '$6 > 0' glimmer_3.out > glimmer_3_plus.out $ awk '$6 < 0' glimmer_3.out > glimmer_3_minus.out 4. 다음을실행하여 CDS 서열을추출한다. plus 및 minus strand 각각을실행할때입출력파일을바꾸는것 말고도어떤점이다른지를유의하여실행하라 ( 파랑색으로표시 ). 아미노산으로번역하려면 EMBOSS 의 transeq 를이용한다. glimmer 가출력하는유전자구조는종결코돈까지를표시하므로번역을마치면마지막 아미노산뒤에 * 이추가되며, 서열 ID 뒤에 _1 이부가된다. 이는 vi 에서수정하면된다. $ awk '{print $1, $3, $4, $5' glimmer_3_plus.out while read seqid ORFid start end > do
> seqret -sbegin $start -send $end -stdout -auto contigs.fa:$seqid sed 's/^>.*$/'">$orfid $seqid:$start-$end +"'/' >> orf_plus.fa > done $ awk '{print $1, $3, $5, $4' glimmer_3_minus.out while read seqid ORFid start end > do > seqret -sbegin $start -send $end -sreverse -stdout -auto contigs.fa:$seqid sed 's/^>.*$/'">$orfid $seqid:$start-$end -"'/' >> orf_minus.fa > done $ cat orf_plus.fa orf_minus.fa > orf.fa $ transeq orf.fa orf_translated.faa 5. seqret 대신 samtools faidx 를이용해서도 EMBOSS 의 transeq 를이용하여같은작업을할수있을것이다. Multiple fasta 파일을서열별로분리하여별도의파일로저장하기 (fastasplit.pl)!/bin/env perl Usage : fastasplit.pl multiple_fasta_file [output_directory] use warnings; use strict; my $inputfile = shift; my $outdir = shift "."; if (! -w $outdir ) { die "You cannot make files in this directory $outdir: $!"; open INPUT, $inputfile or die "Can't open $inputfile: $!"; while ( <INPUT> ) { if ( /^>([^\s^\t]+)/ ) { my $outfile = $1; my $outfilewithpath = "$outdir/$outfile"; open OUT, ">$outfilewithpath"; print STDERR "Writing file $outfilewithpath\n"; print OUT $_; else { print OUT $_; close OUT; Contig 서열파일로부터 assembly 관련수치구하기 (n50.pl)!/usr/bin/perl Source: http://genomics-array.blogspot.com/2011/02/calculating-n50-of-contig-assembly-file.html Modified by Haeyoung Jeong Read Fasta File and compute length my $length; my $totallength; my @arr; my $num;
while(<>){ chomp; if(/^>([^\s]+)/){ print " $length\n" unless $num == 0; $num++; print "SEQ: $1 "; push (@arr, $length); $totallength += $length; $length=0; else { s/\s//g; s/\t//g; $length += length($_); push (@arr, $length); for the last contig print " $length\n"; $totallength += $length; close(fh); my @sort = sort {$b <=> $a @arr; my $n50; print "$num contigs (total length: $totallength)\n"; print "Max. contig length is ", $sort[0], "\n"; foreach my $val(@sort){ $n50+=$val; if($n50 >= $totallength/2){ print "N50 is $val (N50 contig length is $n50)\n"; last; print "Average contig length is ", $totallength / $num, "\n";
참고사이트 [Srrenu Vattipally] Bioinformatics I/O: Tips & tricks from a cluster of bioinformaticians http://bioinformatics.cvr.ac.uk/blog/category/awk/ [Frederick J. Tan] sed, AWK, and bash scripting https://emb.carnegiescience.edu/sites/default/files/140602-sedawkbash.key_.pdf Bioinformatics one-liners https://github.com/stephenturner/oneliners seqtk: Toolkit for processing sequences in FAST/Q formats https://github.com/lh3/seqtk sed and awk for genomics http://macmanes.com/sed-awk/