2022년 3월 30일 수요일

파이썬의 딕셔너리와 펄의 해시를 비교해 보다

펄(Perl)의 데이터 유형 중 하나인 레퍼런스는 어떤 데이터 유형이든 가리킬 수 있다. 그런데 쓰임새가 약간 독특해서 '내가 레퍼런스를 쓰고 있다'는 것을 염두에 두고 코딩을 하지 않으면 혼동을 하기 쉽다. 간혹 anonymous 형태의 자료를 쓰기도 한다.

파이썬 초급 공부를 하면서 딕셔너리(dictionary)라는 데이터 유형은 펄의 해시와 아주 비슷하다는 것을 알게 되었다. 그런데 값을 넣어야 하는 자리에 태연하게 리스트를 넣는 것을 보고 약간의 문화적 충격(?)을 받았다. 간단한 파이썬 예제 코드를 살펴보자.

>>> var = {'kor':[100, 80, 90], 'eng':[88, 65, 44], 'mat':[80, 85, 99]}
>>> type(var)
<class 'dict'>
>>> var['kor']
[100, 80, 90]
>>> var['kor'][0]
100

다음은 펄 코드.

my @array = ( 80, 85, 99 );
my %var = (
    kor => [ 100, 80, 90 ],
    eng => [ 88, 65, 44 ],
    mat => \@array
);
print $var{'kor'}->[1], "\n";
print $var{'mat'}->[2], "\n";

파이썬에서 리스트를 선언할 때에는 [..] 안에 원소를 넣지만 펄은 (..)이다. 그런데 펄에서 배열에 대한 anonymous array를 만들어 reference에 할당할 때에는 [..]를, anonymous hash를 만들어 같은 일을 할 때에는 {..}를 써야 한다.  아, 헷갈린다. 

파이썬의 pandas 라이브러리를 더듬더듬 쓰기 시작하면서 또 걱정 거리가 생겼다. R의 데이터프레임 사용법과 혼동하기 딱 좋으니 말이다. 가끔 용법이 비슷한 주제에 대해서 각 프로그래밍 언어를 비교하는 표를 만들어 암기를 해야 될 것만 같다. 비교 프로그래밍 언어학이 정말로 필요한 순간이다.

2022년 3월 27일 일요일

[왕초보 파이썬 - pandas] 딕셔너리로부터 데이터프레임 만들기

구글에서 조금만 검색을 해 보면 파이썬 초보자를 위한 학습 자료가 넘쳐난다. 더 이상 늦을 수는 없다는 비장한 각오로 파이썬 공부를 시작하면서 왜 블로그에 글을 남기는가? 이것은 '나를 위한 메모'라는 성격이 강한 글이다. 일기는 일기장에 쓰는 것이 맞지만, 여기는 나의 개인 블로그이므로 나의 기억을 돕기 위해서 기록하는 것이니 별 문제가 없을 것이다.

numpy ndarray를 만들어서 데이터프레임을 만들 때에는 특별히 지정하지 않는 경우 행(row) 단위로 결합이 이루어지는 것 같다. pandas.DataFrame 공식 문서(링크)에서 'Construct DataFrame from numpy ndarray' 항목을 찾아보라.

이번에는 4명의 성적 자료를 담은 간단한 데이터프레임을 만들어 보자. 과목은 3가지이다. 각 과목에 대해 개인별 시험 점수를 담은 딕셔너리를 만들고, 이를 pandas의 DataFrame 함수로 합치려 한다. 열 단위로 합칠 것인가, 혹은 행 단위로 합칠 것인가? 기본 동작은 ndarray를 합칠 때와는 달리 열 단위로 합쳐진다. 행 단위로 합치려면 pandas.DataFrame.from_dict(... orient='index'...)를 사용해야 한다.

>>> import pandas as pd
>>> from pandas import DataFrame
>>> data = {'eng':[10, 30, 50, 70], 'kor':[20, 40, 60, 80], 'math':[90, 50, 20, 70]}
>>> names = ['Kim', 'Park', 'Lee', 'Choi']
>>> df = pd.DataFrame(data, index=names)
>>> df
      eng  kor  math
Kim    10   20    90
Park   30   40    50
Lee    50   60    20
Choi   70   80    70
>>> df_2 = pd.DataFrame.from_dict(data, orient='index', columns=names)
>>> df_2
      Kim  Park  Lee  Choi
eng    10    30   50    70
kor    20    40   60    80
math   90    50   20    70

속성(attribute)과 메쏘드를 구별하는 것도 아직 어려운데 데이터 파일을 어느 정도 주무를 수 있게 익히려면 아직 갈 길이 멀었다. 머릿속에서 R 문법과 뒤죽박죽이 되면 어쩐다?


나는 왜 파이썬을 배우려 하는가?

데이터를 분석하고 시각화하기 위해 R을 어느 정도 수준(?)으로 쓰고는 있다. R이 data wrangle, visualization 및 reporting에는 정말 도움이 되는 도구라는 데에는 아무도 이견이 없을 것이다. 하지만 점차 필수가 되어 가는 기계 학습을 공부하기 위해서는 파이썬 환경에서 scikit learn이나 tensorflow를 다루어야만 한다. 또한 생명정보학 분석용으로 공개된 수많은 파이썬 애플리케이션을 제대로 쓰려면 단지 '설치와 활용'을 잘 하는 것으로는 부족하다. 그래서 파이썬이 필요하다.

다행히도 파이썬과 R은 궁합이 잘 맞는다. 다음과 같은 글이 공연히 희망을 갖게 한다. 


"The bottom line is that knowing both R and Python makes you SUPER PRODUCTIVE".

2022년 3월 24일 목요일

구태에서 벗어나기 - Legacy BLAST는 이제 그만!

일루미나 플랫폼으로 라이브러리를 만들어서 시퀀싱을 할 때에는 PhiX174 박테리오파지 DNA가 control로 들어가게 되는데, 이것이 최종 결과물까지 남아서 NCBI에 유전체 정보를 등록할 때에도 그대로 따라가는 경우가 종종 있다. 논문으로도 나올만큼 잘 알려진 문제이다.

Large-scale contamination of microbial isolate genomes by Illumina PhiX control. Standards in Genomic Sciences volume 10, Article number: 18 (2015)

ZGA 파이프라인에는 유전체 조립 결과물에서 PhiX174 박테리오파지 염기서열을 점검하는 기능이 있다. 평소에는 활성화되지 않으므로 다음과 같이 옵션을 주어야 한다.

--check-phix          Check genome for presence of PhiX control sequence.

이 기능을 시험하다가 실제 명령어가 어떻게 이루어졌는지 궁금하여 zga.log 파일을 열어 보았다. 단순하에 blastn을 이용하고 있었다.

2022-03-24 10:15:29,851 - DEBUG - Running: blastn -query /opt/miniconda3/envs/zga/lib/python3.7/site-packages/zga/data/phiX174.fasta -subject /data/project/AA37_KRIBB_Edu_2022-Mar/02_data/_check_phiX/assembly/assembly.fasta -outfmt 6 sseqid pident slen length -evalue 1e-6

그런데 옵션 중 하나가 눈길을 끌었다. '-db'가 아니고 '-subject'라고? makeblastdb로 blast DB를 만들지 않고 단지 FASTA file을 제공해도 된단 말인가? Legacy blast를 즐겨 사용하던 나는 충격을 받았다. 왜 나는 지금까지 legacy blast를 쓰고 있었나? 단지 익숙하기 때문이다. 최신 blast, 즉 blast+는 아직도 옵션을 다 외우질 못해서 매번 도움말을 참조해야 한다. 명령행에서 참조를 하거나 아니면 NCBI 웹사이트를 통해서 말이다.그런데 데이터베이스를 만들지 않고도 blast를 실행할 수 있었단 말인가? 물론 legacy blast에서 bl2seq를 통해서 서열 두 개를 비교할 때에는 DB를 쓰지 않지만, 파일 두 개가 아니라 서열 두 개라는 점이 중요하다.

Blast DB를 만들어서 사용하는 이점이 당연히 있을 것이다. 검색 속도가 빠를 것임은 자명하다. 미생물 유전체 32개를 모아서 하나의 파일(all.fna)을 만든 뒤, 두 가지 방법으로 전부 blastn 검색을 해 보았다.

$ time blastn -query query -subject all.fna -out 1.out

real	0m2.959s
user	0m2.734s
sys	0m0.225s
$ makeblastdb -in all.fna -dbtype nucl -out DB -title custom_DB

Building a new DB, current time: 03/24/2022 11:13:37
New DB name:   /data/project/AA37_KRIBB_Edu_2022-Mar/02_data/09_ANI/fasta/DB
New DB title:  custom_DB
Sequence type: Nucleotide
Deleted existing Nucleotide BLAST database named /data/project/AA37_KRIBB_Edu_2022-Mar/02_data/09_ANI/fasta/DB
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 812 sequences in 0.993721 seconds.

$ time blastn -query query -db DB -out 2.out

real	0m0.192s
user	0m0.149s
sys	0m0.043s

DB를 만드는 쪽이 훨씬 빠르다. 너무나 당연한 이야기겠지만... 결과 파일에서는 database를 표시하는 곳에서 다음과 같은 차이가 있다. 먼저 파일 자체를 이용한 경우(-subject all.fna)는 이러하다.

Database: User specified sequence set (Input: all.fna).
           812 sequences; 131,110,653 total letters

DB를 만들어서 검색을 한 경우는 이렇다.

Database: custom_DB
           812 sequences; 131,110,653 total letters

makeblastdb를 실행할 때 '-title' 옵션의 값으로 제공한 DB의 제목이 여기에 표시된다. 위의 사례에서는 custom_DB였다. 검색 단계에서 DB 타이틀 명칭을 'blastn -db custom_DB'처럼 입력해서는 안 된다. 여기에 들어갈 값은 'formatdb -out' 명령에서 지정한 DB의 이름이다. DB title은 blast 검색 결과 파일에나 표시될 것이다. DB의 타이틀과 이름을 별도로 지정하는 것이 헷갈렸던 것도 최신 blast를 적극적으로 쓰지 않은 구차한 이유 중 하나다. 오늘부터는 그럴 일이 없을 것이다.

하나의 자료에 두 개의 accession, 그러한 것은 있을 수가 없다!

3월 22일에 올린 글(링크)에서 하나의 bacterial genome record에 두 개의 accession number가 붙게 된 이야기를 했었다. 동일 자료에 대해서 CP093535CP093294 두 개의 번호가 붙었다. NCBI 검색창에서 두 번호는 똑같은 수준의 취급을 받지는 않는다. 굵게 표시한 앞 번호만 검색 결과를 내보인다. 그러나 웹브라우저의 주소창에 다음과 같이 입력하면 둘 다 같은 곳으로 연결된다. 그러나 후자의 경우 표시되는 주소가 끝부분이 GI로 치환된다.

  • https://www.ncbi.nlm.nih.gov/nuccore/CP093535
  • https://www.ncbi.nlm.nih.gov/nuccore/CP093294 -> 주소창에 적었던 정보가  https://www.ncbi.nlm.nih.gov/nuccore/2209443742로 바뀜

이 숫자를 GI라고 생각하게 된 것은 바로 다음의 검색 결과 때문이었다.

GI 번호로 검색은 되지만 정작 자료를 열어보면 그 안에는 GI 번호가 없다.

GI는 2016년에 공식적으로 퇴출된 것이 아니었나? 제출 자료가 점차 늘어나고 동일 샘플에 대한 유전체 정보 업데이트도 이어지면서 NCBI는 정수로 이루어진 기존의 GI를 accession.version 형식으로 바꾸기로 하였다. 따라서 flat file이나 FASTA format 등 많은 곳에서 더 이상 GI를 보기 어렵게 되었다.

NCBI is phasing out sequence GIs - Here's what you need to know

2022년 3월 18일에 공개된 유전체 정보를 GI 번호로 검색할 수 있다는 것은 이례적인 일이다. 내부적으로만 쓰고 사용자에게 드러나는 영역에서만 가리워진다는 뜻인지?

CP093294는 웹 접속을 통한 검색에서는 약간의 '차별'을 받지만 명령어 환경에서는 어떤지 알아보았다.

$ ncbi-acc-download CP093535 # CP093535.gbk 정상 다운로드
$ ncbi-acc-download CP093294
Failed to download file with id CP093294 from NCBI
NCBI Entrez returned error code 400, are ID(s) CP093294 valid?
$ esearch -db nuccore -query "CP093535 [ACCN]"
<ENTREZ_DIRECT>
  <Db>nuccore</Db>
  <WebEnv>MCID_623bb246862c0533902fd1df</WebEnv>
  <QueryKey>1</QueryKey>
  <Count>1</Count>
  <Step>1</Step>
$ esearch -db nuccore -query "CP093294 [ACCN]"
<ENTREZ_DIRECT>
  <Db>nuccore</Db>
  <WebEnv>MCID_623bb24f56cb4a102f527a0e</WebEnv>
  <QueryKey>1</QueryKey>
  <Count>0</Count>
  <Step>1</Step>
</ENTREZ_DIRECT>
$ esearch -db nuccore -query "2209443742 [GI]"
<ENTREZ_DIRECT>
  <Db>nuccore</Db>
  <WebEnv>MCID_623bb284b5a53d59fc628304</WebEnv>
  <QueryKey>1</QueryKey>
  <Count>1</Count>
  <Step>1</Step>
</ENTREZ_DIRECT>
$ esearch -db nuccore -query "2209443742 [GI]" | efetch -format gb > download.gbk

CP093294를 이용하여 검색을 하면 명령행 환경에서도 결과는 나오지 않는다. 그러나 GI 번호로는 검색이 잘 되고, GenBank flat file도 정상적으로 다운로드된다. 물론 내용물은 원본에 해당하는 CP093535이다. 흥미로운 것은 GI 번호를 이용하여 반환한 결과물 어디를 보아도 GI 번호 자체는 보이지 않는다. 검색어로는 자격이 있는데 실체는 없는, 말하자면 유령인가? 

EDirect utility로 다운로드한 GenBank flat file. Accession 항목에는 여전히 CP093294가 보인다. 이것이 진짜 accession number인가? 결론적으로는 단순한 '별명'으로 보아야 한다.

Riley 박사가 이미 발급받은 accession number를 제한적으로나마 활용할 수 있게 배려해 주었지만, 결국 하나의 genome record가 두 개의 accession number를 갖게 된 것은 전혀 아니었다. 데이터베이스를 다루는 사람의 입장에서는 이는 허용해서는 안 될 일이다. 공식적이고 유일한 accession에 대한 별명이 하나 붙었다는 정도로 이해하면 된다.

GI는 외부 사용자에게 보이지 않을 뿐 없어진 것은 아니었다. 2209443742을 하나씩 증가시키거나 줄이면서 NCBI에서 검색을 해 보면 감이 잡힐 것이다. 

2022년 3월 22일 화요일

[GenBank] 하나의 염기서열 자료에 두 개의 accession number가 붙었다

작물에 풋마름병(bacterial wilt)을 일으키는 식물병원체인 Ralstonia pseudosolanaceraum 균주의 유전체 정보를 등록하다가 경험한 일을 기록하고자 한다.

신젠타 웹사이트 화면 갈무리(링크).

PacBio 시퀀싱 raw data를 Unicycler로 조립하여 얻은 총 4개의 complete sequence를 GenBank에 올려서 accession number까지 받은 상태에서, 같은 자료를 과거에 HGAP으로 조립하고 polishing을 추가로 실시했던 염기서열의 quality가 훨씬 좋은 것을 발견하였다. 더 낫다는 근거는 PGAP으로 주석화한 뒤 pseudogene의 숫자가 압도적으로 적음에 의한다. 

Accession number는 받았으나 아직 공개가 되기 전이니 염기서열을 바꾸어도 되겠느냐고 이메일을 보냈다. 늘 친절하게 답장을 보내주는 Leigh A. Riley 박사는 그래도 괜찮으니 혹시 예전 accession number를 그대로 쓰고 싶으면 알려달라고 하였다. 마침 이미 확보한 accession number를 언급한 논문을 막 투고한 상태라서 그렇게 해 달라고 하였다.

어떤 학술지는 논문에 기재한 accession number가 투고 당시에 접속 가능한 상태가 되도록 요청하는 경우도 있고, 또 어떤 학술지는 이를 사전에 꼭 요청하지는 않는다. NCBI의 정책은 염기서열 등록자가 요청한 공개 개시일(등록일 기준 최대 3년이었던가?)이 되거나 혹은 그 accession number를 언급한 논문이 출판되는 시점 중 빨리 다가오는 날짜에 맞추어서 공개하는 는 것이다.

그런데 여기에서 문제가 발생하였다. 새로 제출한 염기서열에 새로운 번호가 붙어 버렸고, 이미 확보된 번호로 바꾸는 것이 곤란하다는 것이다. 시차를 두고 이메일 교환을 하다 보면 일어날 수 있는 일이다. 그래도 NCBI의 이메일 응대는 정확하고 빠름을 부정할 수 없다.

대신 Riley 박사는 accession 항목에 기존 것을 병기하는 묘안(?)을 제시하였다. 이렇게 하여 하나의 레코드에 두 개의 accession number가 붙는 기이한 상황이 벌어졌다. GenBank와 RefSeq 전부 이렇게 등록되었다.

GenBank CP093535.1

이런 전례가 과거에도 있었는지는 전혀 모르겠다. 뒤쪽에 붙은 accession number인 CP093294는 NCBI 웹사이트 검색창에 아무리 넣어도 결과가 나오지 않는다. 이상하다? 엊그제까지는 검색이 되었었는데 말이다. 대신 웹브라우저 주소창에 'https://www.ncbi.nlm.nih.gov/nuccore/CP093294'를 입력하면 CP093535의 자료가 잘 나타난다. '.1' 버전번호를 넣으면 안 된다. 입력한 주소는 'https://www.ncbi.nlm.nih.gov/nuccore/2209443742'로 바뀌어 나타난다. 

NCBI 측에서도 나름대로 최선을 다 한 것으로 볼 수 있다. 오늘 아침, 학술지 측에서 accession number에 대한 하이퍼링크를 달아서 원고를 다시 제출해 달라는 연락이 왔다. 기와 이렇게 되었으면 URL이나 검색을 통해서 정상적인 결과가 나오는 나중의 accession 번호를 쓰는 것이 낫겠다 싶어서 수정한 원고를 다시 제출하였고, 곧바로 학술지 측에서 OK라는 답장이 왔다. Entrez Direct로 이 레코드에 대한 accession을 요청하면 어떤 답이 나올지 궁금하다.

참으로 기묘한 경험이었다.

2022년 3월 17일 목요일

sort -V의 힘 - natural sort of (version) numbers within the text

sort는 매우 유용한 명령어이지만 다루기에 까다로운 점도 없지 않다. 버그가 있다는 소문도 있다. GFF3sort라는 Perl 스크립트의 GitHub 사이트에 이런 글이 있다.

However, either the GNU sort or the gt tool has a bug: Lines with the same chromosomes and start positions would be placed randomly. Therefore, parent feature lines might sometimes be placed after their children lines. 

여러 컬럼으로 이루어진 텍스트 파일을 원하는 컬럼 기준으로 정렬할 때 sort는 힘을 발휘한다. 만약 3번째 컬럼을 기준으로 정렬하려면 'sort -k3,3 <filename>'이라고 입력하는 것이 정석이다. '-k 3' 하나만 입력하면 될 것 같은데, 여기에는 약간 복잡한 내용이 숨어 있다. 이 글에서는 그것까지는 상세하게 다루지 않겠다.

1, 2, 3... 10, 11, 12...를 정렬하면 컴퓨터는 1, 11, 12... 이런 식으로 결과를 낸다. 인간이 자연스럽게 느끼는 숫자의 오름차순으로 만들려면 sort 명령어에서 다른 옵션을 주어야 한다. 그런데 정렬 기준이 되는 컬럼이 문자와 숫자의 복합체라고 가정해 보자. 예를 들어서 vsearch로 dereplication한 염기서열 파일의 ID를 생각해 보자.

NR_1;size=2

빨간색으로 표시된 숫자를 기준으로 sort를 하고 싶다. 아예 Perl 스크립트를 짜서 해결할 수도 있지만 sort 명령어 안에 해결책이 숨어 있을 것만 같았다. 구글을 뒤져 보니 StackOverflow에 이런 글이 있었다.

How to sort a file, based on its numerical values for a field?

그것은 바로 'sort -V'를 쓰는 것이었다. 'man sort'에서는 이에 대한 설명이 나오지 않는다. 'info sort'를 해야 겨우 알 수가 있는 기능이다. 이를 사용하니 문자와 숫자가 복잡하게 얽힌 필드를 기준으로 숫자 오름차순의 정렬을 할 수 있었다.

sort, join 등 문자를 직접적으로 다루는 명령어들은 문자열의 사전적 배열 순서 또는 Locale 등에 민감하므로, 그 용법을 정확히 알고 써야 한다. 

2022년 3월 16일 수요일

Contig circularity를 점검하는 스크립트, check_circularity.pl

박테리아 유전체를 long read 기술로 시퀀싱하여 최신 assembler로 조립하면 원형 구조의 replicon 구조를 쉽게 얻는다. Contig가 원형으로 닫힌 구조인지 어떻게 확인하는가? Contig 앞뒤에 나오는 서열이 겹치면 그렇게 판단할 수 있다. 자잘한 contig를 정리하고, 부족한 부분은 채워서 circular contig로 만들고, 마지막으로 dnaA 유전자를 찾아서 염색체 시작 위치를 조정해 주는('fixstart') 종합 도구로는 circlator만한 것이 없다. 아쉽지만 이 프로그램이 앞으로 업데이트될 가능성은 별로 없다고 보는 것이 옳다. roary도 마찬가지일 것이다. 두 프로그램은 모두 2015년에 영국의 Sanger Institute에서 개발하였다.

요즘 나오는 long read용 조립 프로그램(flye나 unicycler가 대중적임)은 대부분 자체적으로 circularization을 해 준다. 단, canu는 예외이다.

Circlator와 같은 종합 도구를 실행하지는 않고 단지 주어진 염기서열의 시작과 끝에 같은 방향으로 반복하여 나타나는 서열이 있는지를 점검하고 싶다고 가정하자. gepard를 써서 시각적으로 확인할 수도 있고, contig 앞부분을 조금 잘라서 전체 염기서열에 대해 blastn을 해도 되지만 너무 귀찮다. 좀 편하게 확인할 길이 없을까?

Circlator 논문을 읽다가 sprai assembler 패키지의 check_circularity.pl이라는 스크립트를 개선하여 additional file 3에 추가하였다는 내용을 보았다. 만약 작은 플라스미드처럼 길이가 read length보다 짧은 경우 tandem repeat 형태의 contig가 나올 수 있는데, 오리지널 check_circularity.pl 스크립트는 이를 1 반복 단위로 줄이지는 못한다. Circlator 패키지에는 당연히 이 기능이 들어 있으며, 친절한 개발자들은 이미 알려진 다른 프로그램의 부속 스크립트도 손을 보아서 논문에 추가해 놓았다. 논문의 'Optimizing existing method' 섹션 일부를 인용해 본다.

출처: Genome Biology(2015)

SPRAI(Single Pass Read Accuracy Improver)는 long reads 시퀀싱 기법이 나온지 얼마 되지 않던 시절 그 이름을 들어 본 일이 있다. Sprai는 지금은 유지되지 않는 웹사이트에만 프로그램이 공개되었고 그 소프트웨어 자체를 발표한 논문은 나온 일이 없다. 나는 sprai가 제공하는 check_circularity.pl 스크립트의 단순함이 마음에 들어서 이를 구해 보기로 하였다. 검색을 거듭한 결과 anaconda package 웹사이트에서 *.bz2 파일의 링크를 겨우 구할 수 있었다. 파이썬 2.7을 쓰는 오래 된 패키지이니 conda environment를 만들어서 설치할 필요까지는 없고, 압축 파일을 받아서 푼 뒤 내가 원하는 스크립트만을 찾아서 shebang line을 수정하였다.

자, 그러면 시작과 끝 부분에 중복된 서열이 나오는 contig를 만들어서 시험해 보아야 한다. 앞서도 이야기했지만 요즘 인기 있는 assembler는 circularization을 다 해서 제공을 하니 check_circularity.pl을 이용할 여지가 없다. 그러나 canu는 다르다. Canu 공식 웹문서에서 제공하는 대장균의 25x PacBio read를 가져다가 canu v2.2로 조립을 하였다.

>tig00000001 len=4665109 reads=10220 class=contig suggestRepeat=no suggestBubble=no suggestCircular=yes trim=20235-4660825

참으로 친절한 설명이다. 이제 check_circularity.pl을 실행해 본다.

$ check_circularity.pl ecoli.contigs.fasta tempdir
tig00000001	circular	[1-24522] matches [4640588-4665109] with 99.833% identity.

서열의 match 여부는 blastn으로 점검한다. 그런데 canu contig의 서열 ID를 열어보다가 예전에는 없었던 'trim=20235-4660825'이라는 정보가 눈에 확 들어온다. 여기에서 제시하는 coordinate만 취하면 서열 앞과 뒤에 겹쳐서 나타나는 영역을 제거할 수 있음이 아니겠는가? 허 참, 세상은 넓고 새로운 것은 없구나. 새로워 보이는 것은 이미 존재하는 것을 뒤늦게 발견한 것에 불과할 뿐...

2022년 3월 14일 월요일

[Perl programming] 배열(array)과 목록(list)은 다르다!

Perl에서 배열(array)과 목록(list)은 종종 같은 것으로 취급된다.  그러나 아주 엄밀하게 말하자면 두 가지는 다른다. Perl에서 배열은 데이터 구조의 일종이다. 즉 스케일러(스칼라? scalar), 배열(array), 그리고 해시(hash)라는 세 가지 데이터 구조 중 하나를 의미한다.

목록은 말 그대로 목록이지 데이터 구조가 아니다. 이는 단순히 스케일러들의 '나열'된 형태를 의미한다. 따라서 전체 목록에 대하여 어떤 변수를 할당한 뒤 인덱스나 키 값에 의하여 그중의 일부를 접근하거나 조작할 수는 없다. 즉 $var = (10, 20, 30);과 같이 쓸 수는 없다.

[GeeksforGeeks] Perl - Arrays vs Lists 

리스트에 변수를 할당할 수는 없지만, 그 자체로서 인덱스/루프/레인지 관련 작업은 가능하다. 얼핏 보기에 기괴하지 않은가? 변수로 할당되지 않은 상태에서 이런 일을 할 수 있다는 것은... 위에서 소개한 GeeksforGeeks에 실린 코드를 살펴보자. 어디에도 변수명 같은 것은 없지만, 이 코드는 잘 돌아간다.

#!/usr/bin/perl

print("Accessing element at index 2: ");
print((10, 20, 30, 40, 50)[2]);
print("\n\n");

print("Range function on list\n");
print join(' ', 1..6);
print("\n\n");

Perl 프로그램에서 괄호와 관련한 다른 궁금증을 해결해 보기로 하였다. 바로 'my'에 관한 것이다. my는 함수인가? 'man perlfunc'(Perl builtin function에 대한 매뉴얼)를 실행하면 다음과 같이 my에 대한 설명이 나오는 것으로 보아서 함수는 맞는 것 같다.

그러나 이 문서의 앞부분('Perl Functions by Category')을 잘 살펴보면 perlfunc 매뉴얼은 Perl에서 함수 또는 '함수처럼 보이는 것'(keyword 및 named unary operator)까지 포함하여 설명함을 알 수 있다. my는 named unary operator의 하나이며(참조) 변수의 범위 설정과 관련한 keyword의 하나이기도 하다. 함수인가, 연산자인가? Perl에서는 그 경계가 모호한 것 같다. 'man perlop'를 입력하여 나오는 매뉴얼의 'Terms and List Operators (Leftward)'를 살펴보면 머릿속은 더 복잡해진다. 마치 경전을 읽듯이, 머릿속을 비우고 깨끗한 마음으로 몇 번을 읽어야 오묘한 진리를 터득하게 될 것 같다. 냉정한 0과 1의 세계를 다루는 프로그래밍 언어에서 이런 모호함이라니? 그래서 Perl이 점점 설 자리를 잃는지도 모르겠으나 - 이는 Perl의 철학인 'there's more than one way to do it'과 일맥상통하는 면이 있음 - 그래서 더 재미가 있는 것 아니겠는가?

다음의 두 웹문서를 읽어보기를 권한다.

그렇게 쉽게 읽히지는 않는, 좀 난해한 내용을 담고 있다. 차라리 StackOverflow 웹사이트에 실렸던 질문과 답을 살펴보면 Perl 프로그래머가 현실적으로 직면하는 궁금증이 무엇인지 더 쉽게 파악할 수 있다. 


질문자의 궁금증은 이러하다. my ($var) =  blah...라는 코드를 종종 보게 되는데, 왜 하나의 변수를 괄호로 둘러싸는지 이해하기 어렵다는 것이다. my ($var1, $var2, $var3) = blah...가 올바른 표현이 아닌지를 묻는 것이다. 둘 다 문법적으로 잘못된 것이 아니다. 다만 문맥에 따라서 해석이 달라지는 것에 유의해야 한다. 프로그램 작성자의 의도를 제대로 반영하도록 코드를 짜는 것이 중요하다. 가장 간단한 사례, 즉 my 뒤에 하나의 변수가 오는 경우만 알아보자.
  • '= 배열;' my $variable이라 적으면 $variable 변수는 배열의 크기를 갖게 된다. my ($variable)로 적으면  배열의 첫 번째 원소를 갖는다.
  • '= 목록;' my $variable이라 적으면 목록의 마지막 원소를, my ($variable)로 적으면 목록의 첫 번째 원소를 갖는다.
와, 이거 도무지 헷갈려서 외우기 참 어렵다. 어떤 경우이든 my 뒤에 오는 하나의 변수를 괄호로 둘러싸면 우변에 있는 자료의 첫 번째 것만을 할당한다는 것을 알았다. 아래의 코드를 마치 부적처럼 계속 바라보고 있으면 원리를 꿰뚫게 될 것이다.

@array = ('A', 'B', 'C');
my ($x) = @array;
print $x, "\n";
# A
my $y = @array;
print $y, "\n";
# 3
my ($w) = qw/ a b c d /;
print $w, "\n";
# a 
my $z = ('a', 'b', 'c', 'd');
print $z, "\n";
# d

2022년 3월 12일 토요일

설악산 지게꾼 임기종 씨의 삶과 실직

무거운 짐을 지게에 싣고 설악산을 부지런히 오르내리던 임기종 씨의 삶은 언론에서 꽤 여러 차례 소개되었다. 그의 '직업'이 관심을 끄는 이유는 아무리 보아도 보수에 비해서 너무나 고되고 힘든 일이기 때문일 것이다. 자본주의 사회에서 누군가에게 필요한 노동을 제공하고, 그 댓가를 받는 것은 서로의 필요성에 의해 자연적이고 자발적으로 교환 조건을 합의하여 이루어짐이 옳다. 그러나 여러 가지의 이유로 인해 '기울어진 책상'에서 합의 아닌 합의를 맺는 경우가 많을 것이다. 혹은 합의조차 없이 그저 관성에 의해 과거에 해 오던 일을 계속하기도 한다.

과거 언론에서 다룬 임기종 씨의 삶은 미담 비슷하게 포장된 측면이 많았다고 생각한다. 냉장고처럼 일반인이 평지에서 몇 명 달라붙어서 옮기기에도 힘든 물건을 지게에 싣고 지팡이와 단련된 그의 몸, 그리고 수십 년의 지게꾼 생활에서 몸에 밴 탁월한 균형감각에 의존해서 위험한 산길을 타고 오른다. 무거운 짐을 지고 가다가 잠시 쉬어 갈 때에는 주변의 야생화나 들새를 사랑스러운 눈길로 바라본다(어쩌면 이것은 프로그램 제작 의도 때문에 좋은 그림을 얻기 위해 만들어진 모습일 수도 있음 - 내 생각). 변변치 못한 수준의 품삯을 받으면서도 그는 힘들다거나 돈이 적다고 불평함이 없이 수십 년 동안 산길을 오르내렸다고 한다. 몸과 정신이 불편한 가족을 둔 그는 어찌보면 사회에서 법률로써 마땅히 보호받아야 할 존재여야 한다. 그럼에도 불구하고 임기종 씨는 푼돈을 모아 홀몸노인이나 장애인을 돕는데 써 왔다고 한다. 그 금액이 지금까지 1억원이 넘는다. 도움을 받아야 할 사람이 더 어려운 이를 위해 도움을 주다니? 미담도 이런 미담이 없을 것이다.

조용히 살던 그가 언론에 몇 차례 소개되니 나 자신도 불편한 맘이 드는 것은 사실이었다. 좀 더 나은 대우를 받아야 했던 것이 아닐까? 수행을 위한 암자에 냉장고라니? 높은 산에서 관광객을 맞는 것이 아니라 수행을 위해 머무르는 것이라면 다소 불편한 생활도 기꺼이 감당해야 하지 않을까? 이런 생각을 갖는 것은 나 혼자만이 아니었을 것이다. 급기야는 청와대 국민청원이나 국립공원관리공단 등에 임 씨가 노동력을 착취당하는 것이 아니냐는 지적이 잇따르게 되었고, 결국 그는 일자리를 잃게 되었다고 한다.

사실을 확인해 보자면 임기종 씨는 국립공원관리공단과 직접 근로 계약을 맺은 사람은 아니다. 따라서 공단에 그러한 일이 부당하다고 호소함은 일단 번지수를 바르게 찾은 것은 아닌 셈이다. 그가 고용 계약서 같은 것을 쓰고 암자의 물건 배달 일을 해 왔을 것 같지는 않다. 그의 입장만으로 본다면, 그저 선의로 일을 해 왔고 품삯도 본인이 결정해 왔다고 한다. 왜냐하면 너무 많이 받으면 마음이 편지 않다는 이유 때문이었다.

결론적으로 이번 논란이 일자 그는 일자리를 잃었다. 올해 65세이니 이제 쉴 때가 되기도 했지만, 임기종 씨가 부당한 착취를 당해 온 것이 아니냐는 논란이 일면서 암자 측에서 부담을 느껴서 더 이상 일을 주지 않기로 했다는 것이다.

여기에서 의문이 생긴다. 산꼭대기에 자리잡은 암자는 생활에 필요한 물건을 끊임없이 산밑에서 가져다 주지 않으면 존립이 불가능하다. 그러니 임 씨가 아니더라도 누군가가 계속 물건을 날라 주어야 한다. 그럼 그 일을 앞으로 누가 한단 말인가? 수행 차원에서 승려가 직접? 혹은 다른 지게꾼 후임자가 더 나은 보수를 받고서? 뒷말이 나오지 않을 수준의 높은 보수를 제시한다면 그것은 임 씨에게 계속 일할 기회를 주는 것이 옳다.

불교에서는 수행 정진을 하는 승려를 위해 뒷바라지를 하는 사람이 있는 것으로 알고 있다. 그것이 일반 신도이거나 혹은 사찰에 소속되어 어느 정도의 지위를 갖춘 사람이 하는지는 잘 모르겠지만, 아마 자원 봉사 정도의 수준에서 이루어질 것으로 생각한다. 외주를 주어 해결하던 암자의 짐 나르기를 이제 사찰 내부에서 자원 봉사 차원에서 실시한다? 이것은 그야말로 지구를 거꾸로 돌리는 행위나 다름이 없다. 과거에는 마땅히 해야 할 일로 여겨졌는지 모르지만, 요즘은 교회에서 신도들의 자발적 참여를 통해 청소와 애찬 봉사를 하는 것도 예전과는 다르게 매우 어렵다는 것을 잘 알고 있다.

컴퓨터 앞에 앉아서 약간의 인터넷 조사와 상상만으로 고발의 글을 쓰기는 참 쉽다. 현장에 가서 조금이라도 상황을 알아보기 전에 함부로 말을 해서는 안 된다고 생각한다. 아마 임 씨나 사찰 관계자 모두 이번 일에 대해서 물어본다면 너무 시달려서 더 이상 언급하기 싫다고 손사레를 칠 것이 뻔하다.

참견하기 좋아하는 우리 모두는 어쩌면 임기종 씨가 명예롭게 은퇴할 기회마저 박탈한 것이 아닌가하는 안타까운 마음이 들었다. 이 일과 관련하여 언론의 책임은 어디까지일까? 임기종 씨를 최근에 소개한 방송 프로그램인 『유 키즈 온 더 블럭』이 책임을 져야 한다는 말까지 있었다.

정보, 관심, 그리고 정의감의 발동. 그것이 행동으로 옮겨질 때 어떤 댓가를 치러야 하는지는 아무도 모른다. 그 댓가가 과연 치를 만한 것이었는지 평가하는 것도 쉽지 않다.

2022년 3월 6일 일요일

스웰덤(Swelldom) 헬스 자전거 Super MTS-9800 고치기 - 수리 완료

왜 링고 스타만 자전거에 올라탄 모습을 하고 있는가? CD로 발매된 Beatles의 <Rubber Soul>에서.

각종 운동 기구의 부품을 전문으로 취급하는 곳이 있어서 고장난 부품을 주문하였다.

신해무역 - 런닝머신 헬스자전거 웨이트머신 부품 전문

새로 구입한 부품과 뼈만 남은 부품을 비교해 보았다. 하나가 다 닳아서 없어질 쯤이 되면, 구동 벨트도 교체가 필요할 것이다. 같은 곳에서 헬스자전거 또는 런닝머신용 벨트도 판매를 한다. 우리가 흔히 런닝머신이라고 부르는 운동 기구의 정식 명칭은 트레드밀(treadmill)이다. 헬스 자전거는? Stationary bike이다. 고정식 자전거라고 해야 될까? 요즘 유행하는 스핀 바이크(spin bike)는 보통 모니터가 없고 클럽 등에서 트레이너의 지도에 따라 단체로 몸을 많이 움직이며 타는 자전거를 말한다.

Stationary bike vs. spin bike: How do they differ?

수리를 시작해 보자! 손에 기름때를 묻힐 각오를 하고....

스프링 장력이 꽤 강해서 풀리에 벨트를 걸려면 약간의 요령이 필요하다.

 

벨트를 성공적으로 건 뒤에 좌우 커버를 씌워서 수리를 완료하였다. 이제 바퀴가 아주 잘 구른다. 조만간 진짜 자전거를 정비해서 밖에서 타야 되겠지만, 상황이 허락하지 않는다면 실내 자전거로 만족하자.

오늘 알라딘에서 구입한 중고 CD The Essential Alan Parsons Project를 들으면서 인증 샷을 찍었다. 노라 존스의 첫 번째 앨범인 <Come Away With Me>가 있길래 사고 싶었으나 커버에 금이 가서 포기하였다.

2022년 3월 1일 화요일

비교 프로그래밍 언어학?

'비교 프로그래밍 언어학'(comparative programming linguistics)이라는 학문이 존재할까? 비교 유전체학(comparative genomics)과 비교 언어학이라는 것이 있으니 비교 프로그래밍 언어학이라는 학문 분야도 존재할 수 있을 것이다. 검색을 해 보니 띄어쓰기를 하지 않고 '비교유전체학' 또는 '비교언어학'이라고 한 낱말로 표기함이 옳은 것 같다.

데이터캠퍼스에서 제공하는 빅데이터 콘서트를 이제 3주째 수강하고 있다. 파이썬에 대한 지식이 필요하므로 같은 곳에서 무료로 제공하는 짤막한 파이썬 강좌를 곁들여 공부하게 되었다. 7차시로 구성되어 있어서 프로그래밍 언어 강좌로는 그렇게 길지는 않은데, 지금까지 꽤 여러 차례에 걸쳐서 파이썬 강좌를 들은 경험이 있어서 어렵게 느껴지지는 않는다. 2014년 KOBIC 차세대 생명정보학 교육(인실리코젠)에서 사용했던 교재를 들추어 보면서 두 강좌의 내용을 상호 보완해 본다. 가 서로 어떠한 어쩌면 이번 기회에 파이썬을 제대로 익히지 않으면 앞으로 기회가 영영 없을지도 모른다는 비장한 각오로 공부에 임하고 있다.

빅데이터 콘서트를 주관하는 김원표 대표는 꽤 오래 전에 온라인 통계 교육 강좌를 통해서 처음 알게 되었다. 인공지능 및 빅 데이터라는 시대의 요청에 맞게 회사의 성격을 잘 맞춰 나간 것으로 보인다.

머릿속에는 Perl과 R에 대한 중급(?) 수준의 지식이 늘 맴돌고 있는데 여기에 파이썬까지 더해지니 슬슬 혼동이 일어나기 시작하였다. 데이터의 형을 가리키는 용어가 다르고, 나열형 변수의 원소를 호출하는 인덱싱 기호도 다르고, 조건문을 쓰는 방법도 그러하고... 비교표를 하나 만들어 두고 늘 참조해야 될 것만 같다.

2000년대에 접어들면서 Perl의 인기는 점점 줄어드는 추세이다. Perl에 미래가 없다고는 할 수 없다. 버전 5.3x에서 곧바로 버전 7로 건너뛰면서 새로운 시대를 준비하고는 있다지만, 앞으로의 Perl 사용자는 과거에 Perl을 썼던 사람이라는 슬픈 전망도 있었다.

Per 7 is v5.21 with different settings. (Announcing Perl 7)

파이썬 초급 강좌에서 데이터의 타입과 이를 다루는 방법을 공부하면서 파이썬이 인기를 얻는 것도 당연하다는 생각이 들었다. 공부를 계속 해 나가면 이러한 생각은 더욱 확고해질 것이다. 그렇다고 해도 Perl에 대한 흥미는 계속 유지하고 있을 것이다. 분석 결과를 정리할 때면 어느새 손가락은 '#!/usr/bin/perl'로 시작하는 shebang line을 타이핑하고 있으니 말이다. 아, '#!/usr/bin/env perl'이 더욱 바람직한 모습이지만 말이다.

스웰덤(Swelldom) 헬스 자전거 Super MTS-9800 고치기

집에서 운동을 하려고 헬스 자전거를 사 놓았다가 빨래 건조대가 되었다는 슬픈 전설을 가끔 보게 된다. 우리집의 헬스 자전거는 아마 2010년 이전에 구입한 것으로 기억한다. 그래도 아들과 딸이 번걸아 열심히 타 주어서 본전은 뽑았다고 생각한다. 이 제품은 다나와 검색(결과)에서도 찾아볼 수 있듯이 여전히 잘 팔리고 있다. 

특별히 고장이 날 만한 물건이 아닌데 갑자기 페달이 몹시 뻑뻑해졌다. 뒤로 돌리면 괜찮은데 앞으로 돌리면 거의 돌아가지 않을 정도였다. 텐션 조절을 하는 케이블이 내부에서 끊어져서 항상 최대 부하가 걸렸을 것이라고 여기고 시간이 되면 한번 열어서 점검을 해 볼 마음을 먹고 있었다.

마침 삼일절 휴일을 맞은 오늘, 거실을 대대적으로 정비할 일이 생겨서 정비를 해 보기로 하였다. 본체를 둘러싼 플라스틱 커버를 어떻게 벗겨야 할까? 몇개의 나사를 풀어봐야 커버의 좌우 면만 분리될 분 크랭크에 걸리고 만다. 그렇다면 크랭크를 풀면 될까? 크랭크 볼트 위치의 캡을 벗겼다. 어라? 소켓 렌치가 없으면 풀지를 못하게 생겼다. 자전거를 정비하는데 필요한 거의 모든 공구를 다 갖고 있음에도 불구하고 헬스 자전거를 분해하지 못한다면 말이 되겠는가? 다른 방법이 있을 것이라 생각하고 궁리를 해 보았다. 

일반적인 스패너로는 풀 수가 없다. 본체 커버를 벗긴 뒤 찍은 사진.
페달을 분리한 다음 커버를 비틀고 씨름을 했더니 크랭크를 풀지 않고도 커버를 벗기는 것이 가능했다. 페달을 비롯하여 자전거에는 왼나사로 조이는 곳이 몇 군데 있어서 주의하지 않으면 안 된다. 왼쪽 페달은 왼나사! 이를 잊으면 곤란하다.

도대체 이렇게 뻑뻑하게 된 이유가 무엇인지 알아 보았다. 아!, V-벨트의 장력을 유지하는 롤러가 완전히 마모되어서 회전을 시작하면 벨트가 자꾸 드럼의 반대편으로 빠지는 것이 원인이었다. 자전거 프레임 바닥에 가득하게 떨어진 가루가 바로 롤러의 잔해였다... 

롤러는 벨트의 이탈을 방지하는 날개가 드럼쪽에만 겨우 남아 있었고, 벨트가 닿는 면은 완전히 하얗게 심만 남아 있었다. 그간 자전거가 얼마나 고생을 했을지 미루어 짐작할 수 있다. 이 부품을 빼 보니 그야말로 '뼈만 남은 상태'가 아닌가.



이 부품은 원래 이런 모습이어야 한다. 
헬스장 관리자가 아니라면 이런 부품을 교체할 일도 없을 것이다. 일반 가정에서 V-벨트 텐션 유지용 롤러가 다 닳을 정도로 실내용 헬스 자전거를 달리지는 않을 터이니 말이다. 어쨌거나 이 부품을 사다가 교체를 하면 되겠다고 생각을 했는데 문득 베어링 축 직경이 과연 맞을지 걱정이 되기 시작하였다. 위 이미지에서 잘린 베어링 규격을 다시 찾아보니 '베어링 6000Z'이라고 한다. 도대체 이게 무슨 규격일까?

Deep Groove Ball Bearing (6000Z)

내경은 10mm라고 하니 잘 맞을 것이다.

봄을 맞아서 자전거로 출퇴근을 슬슬 해 볼 생각을 하고 있었는데, 엉뚱하게도 실내용 헬스 자전거를 보수하는 일부터 시작하게 되었다.