16S rRNA gene sequence로는 해상도가 좋은 트리가 나오기 어려울 것이 자명하다. 그래서 택한 도구는 PhyloSift였다. 이 도구는 원래 shotgun metagenomic read로부터 phylogenetic distribution을 알아보기 위하여 개발된 것이지만, 단일 유전체 어셈블리로부터 universal marker(목록; 원래 40개였으나 실제로 쓰이는 것은 37개)를 찾아낼 수 있으므로 미생물 유전체 조립의 완성도를 평가하는데 쓰이기도 한다. 더욱 좋은 것은 alignment까지 된 상태의 서열이 생성되므로 이를 잘 이용하면 tree를 쉽게 그릴 수 있다. PhyloSift가 metagenome sample로부터 만들어내는 tree와는 다른 것이니 주의가 필요하다.
Bacillus_subtilis_168.fa라는 유전체 파일을 이용하여 PhyloSift를 실행했다고 가정하자(search 및 align mode). 그러면 PS_temp/Bacillus_subtilis_168.fa/ 아래에 alignDir 및 blastDir가 생기고, 그 하위에 각 마커에 대한 파일이 쌓인다. marker_summary.txt는 다음과 같이 생겼다.
각 유전체 서열마다 별도로 작성된 이 파일을 병합하여 정리하면 모든 마커가 전부 확인되었는지의 여부를 쉽게 알 수 있을 것이다. 오늘의 고행은 여기서부터 시작되었다. Perl을 사용하면 어떻게 해서든 결과를 얻을 수 있었겠지만, 왠지 R을 쓰고 싶었다. R을 이용하여 코드를 짜 놓으면 재사용성이 높을 것 같았기 때문이었다.
marker_summary.text file의 정리
각 유전체 서열에 대해서 별도의 디렉토리에 나뉜 marker_summary.txt 파일을 한 곳으로 모으되 균주명이 파일의 제목이 되게 바꾸었다. Shell script를 써서 다음과 같이 구현하였다(짧은 코드인데도 공부할 것이 많아서 2019년 3월 8일 설명과 함께 업데이트). 이 스크립트에서는 완벽하지 않다. 털어내야 할 FASTA 파일의 확장자(여기에서는 fasta)를 실제 상황에 맞게 바꾸어야 하고, 출력파일을 저장할 OUTDIR 환경변수를 미리 선언하고 해당 디렉토리도 만들어 두어야 한다. 그리고 for block의 처음에 나오는 <<<도 이해하기가 조금 어렵다. f 변수에는 aaa/bbb/ccc...라는 형태의 값이 저장되는데, 이를 슬래쉬로 구분한 다음 두번째 것인 bbb를 택하겠다는 의미이다. 만약 <<<쓰지 않으면 cut은 f 변수가 가리키는 파일의 내용을 대상으로 작업을 한다. 이 명령행은 FILE=$(echo $f | cut -d'/' -f2)로 바꾸 실행해도 완전히 동일한 결과를 낸다. 16s_reps_bac과 DNGNGWU* 마커에 대한 결과를 같이 출력하고 싶다면 egrep "16s_reps_bac|DNGNGWU" 명령을 써라.
FILES=`find $1 -name marker_summary.txt -type f | grep blastDir`
for f in $FILES
do
FILE=$(cut -d'/' -f2 <<<$f)
FILE=${FILE%.fasta}
# cp $f ${OUTDIR}/${FILE}
grep DNGNGWU $f > ${OUTDIR}/${FILE}
done
find로 찾아낸 파일명(경로 포함)은 다음과 같다. 빨강색 문자열만을 뽑아내는 것이 포인트이다. cut 명령을 사용하여 슬래쉬(/)를 기준으로 분리한 두번째 문자열을 취한 다음 뒷부분의 .fasta를 제거해야 한다.
PS_temp/Paenibacillus_sp._UNCCL52_GCF_000686825.1.fasta/blastDir/marker_summary.txt
오늘의 시도에서는 파일 내부를 grep으로 뒤져서 DNGNGWU...로 시작하는 표준 마커에 대한 결과만을 모았다. match가 없는 마커의 경우 0이 표시되는 것이 아니라 아예 summary file에 존재하지 않는다. 이것이 나중에 R에서 병합 작업을 하는 것을 약간 어렵게 하였다. 만약 모든 마커에 대한 결과가 다 있었다면, 마커의 순서가 동일한지만 확인한 뒤에 cbind()를 하는 것으로 쉽게 끝났을 것이다.
[R] for loop를 이용한 텍스트 파일의 일괄 입력 및 outer join
이러한 용도의 R 패키지가 이미 개발된 것이 있을지도 모르지만 되도록 R의 기본 기능을 충실하에 이용하여 구현하고자 하였다. 먼저 정리된 파일이 모여있는 디렉토리로 이동하여 R을 실행하였다. 다른 곳에서 R을 실행한 뒤 setwd() 함수를 실행하여 파일 위치를 지정하여도 좋다. 파일의 일괄 입력에서는 Import multiple files to R의 예제를 거의 그대로 따랐다.
가 채워지고, 이는 나중에 0으로 채우면 된다. R에서의 join 혹은 merge에 대한 개념은 여기를 참고하였다. 국문 페이지로는 " [R Friend] R 데이터프레임 결합: rbind(), cbind(), merge()"가 매우 유용하다. merge 함수는 두 개의 데이터프레임을 대상으로 한다. 만약 인덱스로 접근 가능한 데이터프레임이 있는데 총 50개가 있다면 이를 어떻게 합치겠는가? 1번과 2번을 merge하여 임시 데이터프레임을 거쳐서 2번에 저장한다(한꺼번에 저장하는 것이 더 간단할지도 모른다). 그 다음 인덱스를 1 증기시키고 2번과 3번을 merge하여 3번에 저장한다. 이를 인덱스=(전체 데이터프레임 수 - 1)이 될 때까지 for loop를 돌린 다음 종료하면 된다.
list.filenames = list.files()다음으로는 데이터 프레임을 join해야 한다. 기준이 되는 컬럼은 markers이다. outer join을 하였으므로 값이 없는 경우에는 NA
# 특정 확장자를 지닌 파일만 읽어들이려면 list.files(pattern="*.txt")
# create an empty list
list.data=list()
length(list.filenames)
# create a loop to read in your data
for (i in 1:length(list.filenames))
{
# 이중 대괄호 안에 인덱스를 넣어야 데이터프레임을 부르는 것이 된다.
list.data[[i]] = read.table(list.filenames[i],sep=",",header=F)
colnames(list.data[[i]]) = c("markers", list.filenames[i])
}
# 다음의 두 명령은 결과가 약간 다르다.
# 궁금하면 class() 함수를 써 보라.
list.data[1]
list.data[[1]]
# add the names of your data to the list
names(list.data) = list.filenames
# name으로 사용된 파일명에서 확장자를 제거하려면 방법을 따로 알아보아라.
# full outer join
for (i in 1:(length(list.filenames)-1))
{
df = merge(x=list.data[[i]],y=list.data[[i+1]],by="markers",all=TRUE)
list.data[[i+1]] = df
}
df = list.data[[length(list.filenames)]]
df[is.na(df)] = 0
rownames(df) = df[,1]
df = df[,-1]
colnames(df) = list.filenames
# 최종 확인
dim(df)
View(df)
 |
| View(df)의 결과 화면(일부). |
> df.2 = data.frame(df[,-1], row.names=df[,1])
> dim(df.2) # 컬럼 수가 하나 줄어들었다.
결과 점검하기
단 하나라도 마커가 발견되지 않은 유전체는 어느 것일까? 또 어떤 마커가 발견되지 않았을까? 마커의 총 수가 37개라고 해서 컬럼의 합을 가지고 판단의 기준을 삼을 수는 없다. 왜냐하면 주어진 마커가 2회 발견되는 경우도 있을 수 있기 때문이다. apply()함수를 사용하면 이를 매우 쉽게 찾아낼 수 있다.
각 row에 대해서 어느 하나라도 0이 존재하는지를 알아보자. 빨강색으로 표시한 마커가 TRUE를 나타내고 있다. 바로 저 row name만을 반환할 수 없을까?
> apply(df.2 == 0, 1, any) # row에 대하여
> apply(df.2 == 0, 2, any) # column에 대하여
위 명령어를 입력한 다음 눈으로 찾을 수는 없는 노릇이다. 정답은 다음과 같다.
> rownames(df.2)[apply(df.2 == 0, 1, any)]
[1] "DNGNGWU00040"
> colnames(df.2)[apply(df.2 == 0, 2, any)]
[1] "Paenibacillus_polymyxa_TD94_GCF_000520775.1"
오늘 하루가 다 지나갔다! 데이터프레임에서 특정 조건을 만족하는 서브셋을 추출하는 것도 중요하지만, 그 조건에 만족하는 column 혹은 row의 이름을 알아내는 것도 중요하다.
댓글 없음:
댓글 쓰기