2018년 11월 4일 일요일

대청호와 신탄진 장 구경

멀리 가지 않아도 가슴이 시원하도록 탁 트인 곳에서 가을 단풍을 즐길만한 곳이 있으니 그곳은 바로 대청호다. 대전 시민에게 계룡산과 더불어서 탁월한 경관을 제공하는 곳 중의 하나이다. 비록 댐이나 방조제가 갖는 부정적 의미에 대해서도 이제 관심을 가져야 할 시대가 되었지만 말이다.

 



전혀 의도하지 않았으나 우리 부부의 얼굴 사이에 하트 모양의 조형물이 자리잡았다. 그런데 가까이 가서 보면 결코 낭만적이지는 않은 문구가 적혀있다. 내용에 대해서는 상상에 맡긴다.
호수 건너편인 북쪽 약간 높은 곳을 바라보면 팔각정 휴게소가 있다. 그쪽은 지대가 높아서 아마 호수를 바라보는 경치는 더 아름다우리라. 시선을 약간 왼쪽으로 옮기면 가파른 산비탈에 절집이 하나 보인다. 이름은 현암사. 기회가 되면 한번 가 보고 싶다. 계단이 많아서 평소에 등산으로 다져진 체질이 아닌 우리 부부는 마음의 준비를 해야 하리라.

마침 나들이를 나섰던 어제는 신탄진 장날이었다. 큰 도롯가에서 파는 단감을 사려다가 아예 본격적으로 장을 보기로 하였다. 유성장(4일/9일)은 꽤 여러 차례 다녀 보았지만 신탄진장(3일/8일)은 처음이었다. 김치가 다 떨어져서 깍두기 거리를 사기로 했다. 장이 서는 골목은 유성장보다는 약간 넓어서 다니기는 수월하였다. 작정을 하고 장을 보러 간 것이 아니라서 발품을 팔면서 어디에서 무엇을 파는지를 파악하려고 애쓰지는 않았다. 유성장은 워낙 많이 다니다보니 골목 입구에서 얼마를 들어가면 모시 송편을 팔고 어느 골목을 돌면 돼지 껍데기 무침과 닭강정이 있는지를 아는 정도는 된다.

잔치국수와 김밥을 곁들여서 약간 늦은 점심을 먹고 길을 나섰다. 무, 쪽파, 생강, 오이, 양파, 양념닭발, 고춧가루, 그리고 홍시를 샀다. 닭발과 돼지껍데기 무침은 집에서 먹은 적이 없는데 딸아이는 그 맛을 알아서 너무나 좋아한다.

광역시에서 열리는 정겨운 장터는 사람들의 사는 내음이 흠씬 풍겨나는 곳이다. 검정 플라스틱 봉지의 사용을 조금만 더 줄일 수 있다면 더욱 좋을 것이다.

독서 기록 - [전환의 시대] 외 네 권


전환의 시대

  • 박노자 지음
선생님이 학생에게 존댓말을 하는 사회, 입시가 사라지고 명문대학도 사라져서 모두가 무상으로 양질의 교육을 받는 사회, 군대가 모병제가 되는 사회, 아프면 무상으로 치료받고 집이 필요하면 공공 영구임대아파트에 입주하는 사회, '여자 같은 남자'라는 말이 남성에게 최고의 칭찬이 되는 사회, 그리고 주말에 바람 쐬라 평양에 다녀오는 사회.

놀 줄 아는 그들의 반격

  • "꾸준한 딴짓으로 기회를 만드는 세대"
  • 오바라 가즈히로 지음|장은주 옮김

나쁜 뉴스의 나라

  • "우리는 왜 뉴스를 믿지 못하게 되었나"
  • 조윤호 지음

생명의 설계도 게놈 편집의 세계

  • "게놈 편집은 우리와 생명의 미래를 어떻게 바꿀 것인가"
  • NHK 게놈 편집 취재반 지음|이형석 옮김
게으른 생물학자는 논문 대신 이 책을 읽음으로써 CRISPR-Cas9 시스템 개발의 역사와 현황 및 의미를 파악하고자 하였다.

스워브(Swerve) ★★★

  • "나를 계속 넓히며 일하는 사람들의 6가지 비밀"
  • 닉 러브그로브 지음|이지연 옮김
전문가들에게 어떤 분야에 대한 질문을 하면 '그건 제 분야가 아닌데요'하면서 한 발 물러서는 경우를 많이 본다. 나 자신도 그런 경향이 없다고는 말하기 어렵다. 마치 공돌리기를 하듯 소관부처가 아니니 책임을 지기 싫다는 뜻도 되겠다. 평생 직업도, 평생 직장도 없는 시대에 한 우물만 파다가는 하나밖에 모르는 바보가 되기 일쑤이다. 말콤 글래드웰이 주장한 '1만 시간의 법칙'은 어떻게 하느냐고? 단순 계산으로도 평생 예닐곱 차례의 전문가가 될 기회가 충분히 있다. 폭넓은 경험과 도전을 하라.

2018년 11월 2일 금요일

43번 오극관(43 power pentode) 싱글 앰프 프로젝트 - [9] 전원 트랜스 추방

진공관 앰프를 만든다고 해서 반도체 부품을 일절 쓰지 말라는 법은 없다. 진공관 앰프의 특성 혹은 본질에 해당하는 부분을 유지하고 있다면 - 즉 출력관과 출력 트랜스 - 나머지는 현대적인 부품으로 구성을 해도 큰 문제가 없으리라는 것이 나의 철학이다. 이에 대해서 반대하는 사람도 있을 것이다. 진공관 앰프 특유의 소리는 초단관에 기인하는 것이라고. 그래서 진공관으로 프리를 구성하고 반도체로 출력단을 만드는 앰프도 종종 보인다. 그러나 Lenard Audio Institute에 다음과 같은 글이 있다(링크).
The heart of a valve amplifier is the output valves and output transformer.   Only the output valves and output transformer contribute to the unique sound character that makes it distinctly different from a solid-state amp.
그래서 나는 이를 믿기로 했다. 중요한 것은 출력단과 출력 트랜스라는 것.

갖고 있던 전원 트랜스를 이용하여 B+ 전원과 히터 전원까지 감당을 하려니 트랜스를 두 개나 써야 하고 덩달아 배전압 정류회로까지 만들어야 하는 불편함이 있었다. 초단/드라이브단에 OP amp 회로를 쓰기로 결정하고 나니 기왕 이렇게 된 것, DC-DC 스텝 업 부스터 모듈을 적극적으로 활용하기로 했다. 그래서 다음 그림과 같은 엽기적인 구성의 앰프가 되고 말았다. SMPS 기반의 어댑터와 부스터가 만들어내는 높은 주파수의 노이즈가 범벅인 상태일지도 모르겠지만 가청 주파수 범위가 아니라면 무슨 문제가 있겠는가?

MT3608을 사용한 부스트 모듈은 바로 어제 국내에서 구입하였다. 주문에서 배송 완료까지 24 시간이 걸리지 않았고, 가격은 하나에 1,100원이었다. 이번 전원부 개조작업의 핵심 역할을 한 고전압 부스트 모듈을 eBay에서 구입하였는데 주문에서 배송까지 무려 45일 가까이 걸렸다. 이렇게까지 긴 시간을 기다리면서까지 취미 제작을 해야 하는가 하는 회의감이 들기도 했다.

LF353 OP amp를 이용한 증폭회로는 원래 +/- 15 V 정도(최대 18 V)의 양전원을 공급해야 한다. 그러나 이는 매우 번거로운 일이므로 단전원 어댑터를 이용하되 저항과 캐패시터로 이를 반으로 나누어서 간이 양전원을 공급하고 있다.

OP amp를 이용한 전압증폭회로에서는 저항의 비율을 이용하여 수백배 이상의 이득을 얻을 수 있다. 하지만 당연하게도 전원전압을 넘어갈 수는 없다. 제이앨범 매니저의 설계에 의하면 43 오극관 싱글 회로를 적정 출력으로 구동하려면 40 V 가까운 범위를 스윙해야 하는데, 내가 마련한 OP amp 회로에서는 출력이 약간 부족할 수밖에 없다. 이에 따라서 제이앨범의 힌트에 따라서 오디오용 트랜스포머를 사용해 본 것이다. 신호의 크기를 가감없이 전달하려면 600 Ω : 600 Ω 트랜스포머를 쓰는 것이 정석이겠으나, 조금이라도 증폭을 하기 위해서 IPT-14라는 것을 뒤집어서 사용하였다. 원리적으로 주파수 특성, 특히 저음부의 전달에서 손실이 있을지는 모르겠으나 귀로 듣기에는 아주 만족스럽다. 출력단을 위한 12 V 어댑터로 OP amp 드라이브단까지 전원을 공급하고, 권선비가 더 높은 트랜스포머를 쓰면 더욱 간단하게 마무리를 할 수 있을지도 모른다.

국내의 웹사이트에서 IPT(input transformer)와 OPT(output transformer)라는 것을 찾으면 두어 가지 종류의 것이 나온다(예: IC114). 그런데 그 용도를 도대체 모르겠다. 전자회로 실습용으로 만드는 전자새 키트 혹은 라디오 키트 말고는 이 부품들이 쓰이는 곳이 보이질 않는다. 권선비와 임피던스비도 잘 이해가 가지 않는다. 이에 대해서는 추후에 좀 더 조사를 해 볼 예정이다.

어찌되었든 오늘은 그림에서 소개한 다이어그램대로 구성을 마쳤고, 무난한 성능을 보였다. 25 V를 만드는 부스터 모듈이 다소 뜨겁기는 한데 55도를 넘을 것 같지는 않다. 며칠동안 테스트를 하여 문제가 발생하지 않는지 확인한 다음 적절한 섀시를 꾸며보도록 하자.



11월 4일에 추가한 글

MT3068 모듈을 잘못 건드려서 두 개를 망가뜨렸다. 12 V 어댑터의 자체 보호 기능이 있어서 망가진 모듈이 단락 상태가 되었지만 자동으로 전원이 차단되어서 더 큰 사고를 예방할 수 있었다. 만약 내가 만든 SMPS를 사용했다면? 아무런 보호 기능이 없어서 어디선가 연기가 나기 시작했을 것이다. 반도체 소자는 무엇인가 잘못되어 파괴되면 도통 상태가 되어버린다는 것이 문제이다.

별도의 24 V 1.5 A어댑터를 사용하여 프리앰프와 진공관의 히터를 연결해 두었다.

용량이 큰 24 V 어댑터가 있다면 히터와 프리앰프에는 직접 전원 공급을 하고, 부스터를 이용하여 고전압을 제공할 수 있을 것이라는 생각을 하지 않았던 것은 아니다. 며칠 전 모니터용 24 V 어댑터를 구해서 이러한 구상대로 연결을 해 보았지만 '웅-'하는 심한 잡음이 발생하였다. 이미 갖고 있던 전원 트랜스와 정류회로를 조합하여 직류를 만들어 보았지만 이것 역시 만족스럽지 않았다. 대용량 캐패시터를 아무리 덧대어 보아도 소용이 없었다.

앰프 전체를 하나의 어댑터로 구동하려면 프리앰프의 출력은 IPT로 아이솔레이션해야 된다. 왜냐하면 OP amp에 공급되는 양전원은 제대로 만들어진 것이 아니라 단전원을 분할한 것이기 때문이다. 따라서 공급전압의 1/2에 해당하는 곳을 그라운드로 삼아서 움직이게 된다. 고전압을 만들어내는 부스터는 입력전압쪽과 분리되어있지 않다. 따라서 진공관 전력증폭회로의 그라운드는 전원의 음극 및 OP amp로 공급되는 직류의 그라운드(신호의 그라운드가 아니라)와 연결된 상태이다. 상황이 이러하니 OP amp 출력의 그라운드를 메인 회로에 연결하는 것은 옳지 않다. 그러나 IPT를 중간에 넣으면 이러한 문제가 해결된다.

12 V 5 A 어댑터와 부스터로 155 V를 만들고 히터(24 V)는 다른 어댑터로 점화하는 것이 현재까지의 최선이다. 전체 구성도는 다음 그림과 같이 바뀌었다.

프리앰프보드를 너무 만지작거리다가 그라운드 신호선의 납땜이 떨어져서 스피커가 부서질듯 '두두두두...'하고 울리는 소음을 접하기도 했다. 출력단자에 거의 전원전압 그대로가 검출되어서 처음에는 OP amp가 망가진 줄로만 알았다. 자작 취미를 하면서 일어날 수 있는 나쁜 모든 일은 모두 다 경험하는 것 같다.

11월 28일에 추가한 글

MT3608 모듈은 한번에 4개를 구입해 놓았었는데 바로 어제까지 전부 망가지고 말았다. 납땜을 하면서 칩에 손상이 올지도 모른다고 생각하여 마지막 것은 볼트로 결선을 했음에도 불구하고 이런 결과를 낳았다. 25V 0.6V를 출력으로 뽑는 것이 그렇게 무리였을까? 분명 specification에서 지정한 수치 이내의 조건이라고 생각한다. 칩과 코일에서는 손을 대기 어려울 정도의 열이 나지만 특별히 방열판을 달거나(달 수 있는 모양이 아니다) 대책을 세우라는 지시는 없었다. 다른 용도로 이 제품을 사용하는 사람들은 어떤지 모르겠지만 최소한 나에게는 좋은 경험이 아니었다.

2018년 10월 27일 토요일

43번 오극관(43 power pentode) 싱글 앰프 프로젝트 - [8] 프리앰프 바꾸기

나무판 위에 서툴게 만든 43 오극관 싱글 엔디드 앰프를 매일 흐뭇하게 바라보면서 어떻게 섀시를 만들지 고민을 하던 중 프리앰프가 망가지고 말았다. 재생 중에 갑자가 전원이 들락날락거리기 시작한 것이다. '퍽' 소리와 함께 파일럿 LED도 꺼지고 아무런 소리가 나지 않는다. 잠시 살아나는 것 같다가 또 꺼지고를 반복하였다. 처음에는 전원 커넥터쪽의 접촉 불량이라고 생각했지만 납땜을 새로 해고 아무리 만져 보아도 소용이 없었다. 전원을 넣고 무슨 일이 벌어지는지 기다리는데 갑자기 연기가 나면서 대용량 저항이 타고 말았다. 문제점을 진단하고 고칠 능력이 없으니 이제 이 앰프는 12AU7 진공관 한 알을 남기고 잡동사니 상자에 들어가는 신세가 되었다.

43 오극관은 전력증폭회로를 담당한다. 이 앞에서 신호(오디오 신호이므로 교류)를 적당한 수준으로 흔들어 줄 '드라이버'가 필요한 것이다. 드라이브단은 소스 기기에서 출력되는 전압을 흔들어 주는 것이라서 큰 전력을 소비하지는 않는다. 명색이 진공관 앰프이니 드라이브단도 진공관을 쓰는 것이 자연스럽다. 그러나 43 오극관의 히터(25 V) 및 B+ 전원(95~160 V)은 요즘 보통 사용하는 다른 진공관과 같이 쓰기에는 적합하지가 않다. 그래서 나는 진공관-MOSFET 하이브리드 헤드폰 앰프 겸 프리앰프를 사용했던 것이다.

제이앨범의 매니저께서는 요즘 43번 오극관을 이용한 앰프 회로의 설계를 거의 마쳐가고 있다. 초단에는 오극관인 6KT6을 사용한다. 계산 결과에 의하면 드라이브단에서는 36~40 V 정도의 범위를 흔들어야 한다. 반도체 프리앰프로는 이는 도달하기 어려운 값이다. 왜냐하면 전원전압이 결국 프리앰프의 출력 범위를 결정하기 때문이다(물론 게인도 중요하다).

이에 나는 수년 전에 만들었던 op amp 기반의 CMoy 헤드폰 앰프(위키피디아)를 개량하기로 하였다. 이것은 헤드폰 앰프 스테이션(하스)의 신정섭 님(sijosae)이 상세한 제작법을 소개하여 국내에 널리 알려진 것이다(제작 정보). 놀랍게도 이 글은 2002년 작성된 글이고 2006년 이후 신정섭 님은 하스 커뮤니티의 자작방에 글을 올리지 않고 있다. 근황이 무척 궁금하지만 나는 하스의 회원이 아니니 물어볼 수가 없다. 그의 작품은 외국에서도 많이 알려져 있다. 궁금하다면 Sijosae's DIY Gallery를 방문해 보자.

헤드폰 앰프는 프리앰프로 써도 별 문제가 없는 것으로 안다(그 역은 항상 성립하지는 않는다). 그렇다면 별도의 목적을 갖고 있는 프리앰프가 진공관 싱글 앰프의 초단에 쓰여도 되는가? 진공관 전력증폭회로가 필요오하는 전압의 폭보다 훨씬 낮은 것은 사실이다. 그러나 침실에서 그다지 크지 않은 음량으로 듣는 것이 목적이라서 문제는 없을 것으로 생각한다.

초단, 드라이브단, 위상반전단...

여기에서 잠시 진공관 앰프의 회로 구성에 대해서 잠깐 논하기로 한다. 위상반전단(inverter)은 푸시-풀 앰프에서 필요한 것이니 제외하기로 한다. 

사진 출처: Lenard Audio Institute
위의 그림을 보면 출력관을 구동(drive)할 전압을 만들기 위해 pre-amp와 driver의 두 단계(stage)를 거친다. 이것은 흔히 이야기하는 초단과 드라이브단을 각각 일컫는 것으로 보인다. 그러나 우리나라에서는 pre-amp라고 하면 특정한 목적이 있는 별도의 장비(즉 power amplifier 혹은 integrated amplifier와 구별되는)로 인식되는 것 같다. 소출력 앰프라면 진공관 한 알을 가지고서 프리앰프과 드라이버의 두 stage를 전부 담당하도록 설계하는 것이 보통이다. 초단, 위상반전단, 드라이브단 및 출력단의 의미를 서병익오디오 기술칼럼에서 확인해 볼 수 있다(링크).

헤드폰 앰프의 개조

내가 신정섭 님의 제작 정보를 따라서 만든 헤드폰 앰프는 게인이 약 1.3 정도에 불과하다. 그래서 저항을 각 채널에서 하나씩 바꾸어서 11배(20.8 dB)로 만들었다. 이렇게 하면 43 오극관 전력증폭회로에 잠시 물려서 사용했던 프리앰프의 그것과 거의 비슷하다. 사용된 op amp는 LF353(datasheet)이다. 


원래 op amp는 직류 양전원이 필요하다. CMoy 헤드폰 앰프에서는 9 V 전지 두 개를 직렬로 연결하여 +9/0/-9 V를 얻는 것이 기본이고, 신정섭 님의 버전에서는 9 V 전지 하나에 캐패시터와 저항을 연결하여 이를 반분하여 +4.5/0/-4.5 V를 만든다. 나는 갖고 있는 SMPS 어댑터를 사용하기로 했다. 9 V, 24 V 및 32 V의 것을 전부 연결해 보았다. 부품의 내압에는 문제가 없다. LF353 데이터시트에 의하면 최대 공급 전압 범위는 ±18 V이니 안전하다. 물론 게인을 더 올리지 않으면 출력 전압이 여기에 접근하지 못할 것이다.

TDA7297 '반찬통 앰프'가 희생되었다.
소리를 들어 보았다. 음량 수준은 바로 수일 전까지 사용하던 진공관-MOSFET 하이브리드 프리앰프와 비슷하다. 고급 op amp라고 할 수는 없지만 잡음 수준도 양호하다.

반찬통 안의 남은 공간이 넉넉하다.

여분으로 보유한 op amp (TL072와 NE5532). 회로 기판에 지금은 LF353이 꽂혀 있는 것은 아마도 이들 중에서 잡음이 제일 적어서 선택된 것 같다.
만약 내가 인내심이 충분했더라면 AliExpress 혹은 eBay에서 몇 달러에 판매하는 NE5532 preamplifier board를 구입하여 사용했을지도 모른다. 그러나 가장 최근에 주문한 물건이 한 달이 훨씬 넘도록 배송이 되고 있지를 않아서 당분간은 해외 구매를 자제하려고 한다.

2018년 10월 26일 금요일

awk와 sed를 사용하여 muti-fasta 파일의 서열 ID를 일괄적으로 바꾸기

수십 개의 대장균 유래 MiSeq 데이터를 Unicycler(GitHub 링크)로 조립하였다. 실행 결과물은 전부 샘플별 디렉토리에 나뉘어 저장된 상태이다. 디렉토리 이름은 전부 다르지만, 그 내용물의 파일 이름은 전부 같다. 예를 들자면 최종 조립 결과 파일의 이름은 모두 assembly.fasta이다. 이를 한데 모아서 저장하자면 파일 이름에 샘플의 이름이 들어가도록 바꾸어야만 한다.

이 작업은 파일의 내용을 직접 건드리는 것이 아니라서 비교적 쉽다.

1
2
3
4
5
$ mkdir final
$ find . -name assembly.fasta | while read f
> do
> cp $f `echo $f | awk -F/ '{sub(/assembly/, ""); print "final/"$2"_plasmid"$3}'`
> done

이상의 명령을 실행하면 ./123/assembly.fasta 파일은 ./final/123_plasmid.fasta로 바뀐다. 세 번째 줄을 참고하면 awk의 sub() 함수를 사용하여 assembly라는 문자열을 제거하고, 파일 경로로부터 샘플의 이름(원래는 디렉토리 이름이었음)을 추출하는 것이 핵심이다.

다음으로는 각 fasta 파일로 들어가서 서열 ID를 바꾸는 방법을 알아보자.

1
2
3
4
5
6
7
$ ls *fasta
1039_plasmid.fasta  1536_plasmid.fasta ... 
$ ls *.fasta | while read f
> do 
> h=${f%%_*}
> sed -i "s/>/>s${h}p/" $f
> done

현 디렉토리에 존재하는 .fasta 파일의 이름을 변수 f에 넣은 다음 밑줄('_')부터 끝까지를 제거하는 명령(다섯 번째 줄)을 먼저 눈여겨 보자. Shell에서 문자열을 조작하는 매우 유용한 방법이니 잘 기억해 두어야 한다. %%, %, ##, #에 따라서 삭제의 범위와 위치 기준이 달라지고, 치환도 가능하다. 상세한 사항은 Advances Bash-scripting Guide: Manipulating Variables를 참조하자. 다음에는 서열 ID를 sed 명령으로 치환한다. sed의 -i 옵션은 파일을 직접 고치는 것을 뜻한다. 따라서 실수를 하면 곤란하다. 그러나 -i.bak라고 입력하면 원본파일.bak이라는 백업본이 남는다.  흥미롭게도 -i와 SUFFIX 문자열 사이에는 공백이 없어야 한다. Shell variable을 sed의 치환 명령어 내에서 그대로 쓸 수 있다는 것이 매우 편리하다. 단, 명령어를 겹 따옴표로 둘러싸야 한다. 이상을 실행하면 서열의 ID는 다음과 같이 바뀐다.

(변경 전) >1 length=25428 depth=1.01x
(변경 후) >s1039p1 length=25428 depth=1.01x

awk를 사용하면 똑같은 일을 좀 더 복잡하게(?)할 수 있다. Shell variable을 awk 명령어 내에서 쓰려면 -v 옵션을 이용해야 하고, 서열 ID가 있는 라인과 서열만 포함한 라인을 모두 출력하려면 next 명령을 써야 한다. 그러나 printf() 함수를 이용하여 좀 더 세련된 결과를 얻을 수 있다. 다음의 사례에서 세번째 awk 명령어에 주목해 보자. awk의 sub( ) 함수 용법을 먼저 이해해야 하고, next의 의미는 곧 이어서 설명하겠다.


1
2
3
4
$ ls *fasta | while read f
> do
> awk -v s="${f%%_*}" '$1~/^>/{sub (/>/, ""); printf ">s%d-p%d %s %s\n", s, $1, $2, $3; next}{print}' ${f} > ./edit/${f}; done
> done

세 번째 줄의 awk 명령어 구조가 좀 난해하다. 이를 이해하려면 다음의 두 명령어가 무엇이 다른지를 알아야 한다.

  1. awk '조건{명령어A}{명령어B}' 파일
  2. awk '조건{명령어A; next}{명령어B}' 파일
작업 중이 라인이 조건에 맞으면 {명령어A} 블록을 실행하고, 맞지 않으면 {명령어B} 블록을 실행하는 것이 아니다! 첫 번째 명령어는 조건에 맞으면 A를 실행하고, 같은 라인에 대해서 그대로 B를 수행함을 뜻한다. B를 수행함에 있어서는 앞에서 선언한 조건과는 상관이 없다. 따라서 조건에 맞는 줄에 대해서 A를 실행하고, 다음 줄로 넘어가려면 그 블록 안에서 next 명령을 선언해야 한다. 이렇게 해야 조건에 맞지 않는 줄에 한하여 명령어 B 블록이 수행된다. 이렇게 표현하는 것이 마음에 들지 않으면 if ~ else ~ 구문을 써야 한다.

아무 fasta 파일에 대해서 다음의 명령을 실행시켜 보라. NR은 현재 작업 중인 라인의 번호를 의미한다. 그 다음에는 next를 제거한 뒤 다시 똑같은 명령을 실행해 보라. 어떤 차이가 있는지를 쉽게 이해할 수 있을 것이다.

$ awk '$1~/^>/{print NR, ": header"; next}{print NR, ": seq"}' test.fasta

오늘 소개한 방법이 최선이 아닐 수도 있다. 다음 웹페이지에서 소개하는 짤막한 awk 코드에도 공부할 것이 많이 숨어 있다.

one liner to split a multifasta into separate single file


1
2
3
4
print 'hello world!'cat hg18.fa | awk '{
        if (substr($0, 1, 1)==">") {filename=(substr($0,2) ".fa")}
        print $0 > filename
}'


2018년 10월 24일 수요일

블로그 포스트에 코드 조각 삽입하기

코드 조각(code snippet)을 블로그 포스트에 예쁘게 삽입하려면 정성을 기울여야 한다. 한때 Syntax Highlighter를 써 보려고 노력했지만 블로그의 템플릿을 수정하는 것도 까다로웠고 실제 적용한 뒤에 페이지를 로드해 보면 제대로 불러오는데 시간이 많이 걸린다는 느낌이 들었다.

좀 더 간편한 방법은 없을까? Stack Overflow에서 그 해답을 찾았다.

Formatting code snippets for blogging on Blogger

첫번째는 hilite.me 웹사이트를 이용하는 것이다. 웹사이트 주소도 아주 외우기 쉽다. '하이라이트 미' 아니겠는가? 여기를 방문하여 소스코드 창에 원하는 코드 조각을 입력하고 언어의 종류와 스타일을 결정한 뒤 나오는 HTML 코드를 복사하여 블로그 편집창에 넣으면 된다. 이렇게 말이다. 편집모드는 HTML로 놓는 것이 좋다.


#!/usr/bin/perl

while (<>) {
    chomp;
    my @data = split /\t/, $_;
    print "$data[0]: $data[3]\n";
}

이러한 용도의 웹사이트는 몇 군데가 더 있다.
또 다른 사용자는 다음의 css 스크립트를 복사하여 쓰는 것을 제안하였다. 중간의 블록에 코드를 삽입한 다음 이를 그대로 블로그 편집창에 넣으라는 것이다. syntax highlighting 기능은 없지만 소스 코드를 그대로 보이는 데에는 유용하다. 매우 단순한 방법이다.

 <pre style="font-family: Andale Mono, Lucida Console, Monaco, fixed, monospace; 
                color: #000000; background-color: #eee;
                font-size: 12px; border: 1px dashed #999999;
                line-height: 14px; padding: 5px; 
                overflow: auto; width: 100%">
       <code style="color:#000000;word-wrap:normal;">

            <<<<<<<YOUR CODE HERE>>>>>>>

       </code>
 </pre>

이를 한 번 적용해 보았다.

#!/usr/bin/perl

while (<>) {
    chomp;
    my @data = split /\t/, $_;
    print "$data[0]: $data[3]\n";
}       

블로그 문서 내에서 '<'과 '>'을 제대로 표시되게 하는 것은 여전히 어렵다. 몇 번씩 '미리보기'로 확인을 해야 하니까 말이다. 오늘 찾은 방법은 잘 기억해 두었다가 활용하도록 하자.

2018년 10월 23일 화요일

Protein sequence clustering과 dN/dS 분석

새로운 과제를 시작하면서 단백질 서열 클러스터링이라는 매우 고전적인 일을 하게 되었다. 단백질 패밀리의 분석, pan-genome 계산, ortholog 추출, dN/dS 분석, 데이터셋의 크기 축소 등 단백질 클러스터링은 매우 여러 분야에 쓰인다. 자료를 조사하면서 지난 2월에 작성했던 글인 Usearch 맛보기도 업데이트하였다. 오늘 작성하는 글은 "How to cluster protein sequences: tools, tips and commands"을 많이 참고하였다.

단백질 서열 클러스터링 기법은 크게 다음의 두 가지로 나눌 수 있다.
  • Sequence-based clustering: CD-HIT, UCLUST, USEARCH...
  • Graph-based clustering: LAST, legacy BLAST, BLAST+...
  • Network clustering: MCL
Graph-based clustering에서는 모든 단백질 서열 사이의 1:1 similarity 정보를 얻어내는 것이 필요하다. BLAST 계열의 프로그램을 사용하는 경우 결과물을 파싱하여 Reciprocal Best BLAST hit(RBB, RBH, or biodirectional best hit BBH)을 추출하여 활용하면 되는 것이다. RBB는 두 genome 사이에서 ortholog candidate를 산출하는 데에도 널리 이용된다.

RBH으로부터 ortholog를 정확하기 추출하기 위한 BLAST option을 제안한 논문이 2008년에 발표된바 있다(PMID: 18042555). 이를 기반으로 하여 다음의 BLASTP 옵션이 제안되었다. 논문에서는 첫번째 조건만을 제안하였다. 이렇게 함으로써 ortholog의 수는 가장 많게, 그러나 error는 최소화할 수 있었다고 한다.
  • -F "m S" -s T (combination of soft filtering and Smith-Waterman final alignment)
  • E value: 10-5 또는 10-6 이하
BLAST 프로그램은 database search와 alignment의 두 단계로 이루어지는데, soft filtering이란 첫번째 단계에서만 low-complexity sequence를 마스킹하는 것을 의미한다(-F "m S"). 기본 조건은 -F T -s F이다. -F "m S"라는 옵션이 있다는 것은 blastall의 모든 옵션이 수록된 페이지(링크; 혹은 명령행에서 blastall이라고만 치면 나오는 텍스트)에서 찾아보기 어렵다.

이상의 옵션은 legacy BLAST에 대한 것이다. 그러면 요즘 권장되는 BLAST+에서는 무엇이라고 옵션을 제공해야 하는가? NCBI의 매뉴얼 페이지(링크)를 참조해 보면 -soft_masking T -use_sw_tback (-evalue 10-5)일 것으로 추청된다. 

다음의 북 챕터에서는 BLAST+를 이용하여 상동유전자 정보를 찾고, 이로부터 dN/dS 분석일 실시하는 방법을 친절하게 설명하였다. 사용된 Perl script(BLAST+ 결과물에서 RBB를 추출하는 스크립트 포함)는 제1저자인 Daniel Jeffares의 웹사이트에서 입수 가능하다고 해 놓았다.

A Beginners Guide to Estimating the Non-synonymous to Synonymous Rate Ratio of all Protein-Coding Genes in a Genome. Methods Mol Biol. 2015;1201:65-90. doi: 10.1007/978-1-4939-1438-8_4. PMID.

그러나 지금 이 사이트를 방문해보면 figshare에 Perl 스크립트 일부만이 남아있다. 그러나 Daniel은 이를 기반으로 하여 VESPA라는 파이프라인을 개발하여 2017년 PeerJ Comput. Sci.에 발표하였다(VESPA: Very large-scale Evolutionary and Selective Pressure Analyses. 저널 GitHub 매뉴얼). VESPA를 공부하려면 단순히 RBB를 찾으려는 목적이 아니라 더 큰 안목을 가지고 착수해야 할 것이다.

InParanoid를 사용해서 두 단백질 세트 사이에서 ortholog를 찾을 수도 있다. 그런데 Why InParanoid sucks라는 글(2014년도)을 보니 이제는 이 프로그램을 떠나야 하는 것이 아닐까하는 생각도 들기 시작하였다. 알고리즘적으로 특별히 문제가 되는 것은 아니지만 말이다. 어쩌면 blast_rbh.py 스크립트가 바로 내가 찾던 것이 아닌지 모르겠다. 이를 다운로드하여 실행해 보았다. duplicate가 있다는 경고는 나오지만 결과 파일에서 이들이 어느 것인지를 보여주지는 않는다. 상동유전자를 찾는 프로그램으로서 GET_HOMOLOGUES도 유용한 프로그램으로 여겨진다.


$ ls
A.faa  B.faa
$ python ~/script/blast_rbh.py -o example.tab -t blastp -a prot A.faa B.faa 

Done, 10 RBH found
Warning: Sequences with tied best hits found, you may have duplicates/clusters
$ ls
A.faa  B.faa  example.tab
$ head -n 1 example.tab 
#A_id B_id A_length B_length A_qcovhsp B_qcovhsp length pident bitscore

BLAST+ 분석 결과물에서 best hit(RBB는 아님)을 추출한 뒤 mcl로 클러스터링하는 명령행을 정리하여 보았다. 이는 서두에서 언급한 참고문헌에서 발췌한 것이다.

$ makeblastdb -in DB.fasta -dbtype prot -out DB
$ blastp -query example.fasta -db DB -outfmt '6 std qlen slen' -out output_tmp.list -evalue 1.03e-05
$ awk '{if($1!~/^#/ && $3>=50) print $1"\t"$2"\t"($3/100)}' output_tmp.list > hits.list
$ sort -nrk 3 hits.list | sort -k 1,2 -u > best_hit.list
$ mcl hits.list -I 1.8 --abc -o clusters.txt
$ awk '{if(NF>=2) print}' clusters.txt > clusters_after_threshold.txt

맨 아래의 awk 명령어는 singleton을 제거하기 위함이다.

마지막으로 대용량의 pan-genome 분석을 위한 새로운 프로그램을 하나 소개해 본다. 늘 Roary만 편식을 심하게 했었는데 이제는 다른 것에도 눈을 돌려볼 때이다.