참고블로그는 다음과 같다.
이 블로그는 위 블로그 내용에 따라 직접 실습한 내용을 적어놓은 학습노트입니다.고주온박사님께 감사드립니다.
자신의 컴퓨터에 설치된 BLAST+를 이용하여 한 번에 분석 결과를 만드는 예를 보자. 초파리 단백질 데이터와 blastp를 이용하여 예제를 진행한다.우선 질의 서열을 몇 개 준비한다. 사용 중인 ~/BLAST/ 디렉터리 밑에 query라는 디렉터리를 만들고, NCBI 웹사이트에서 관심 있는 단백질 서열들을 내려 받아 ~/BLAST/query 디렉터리에 풀어 놓는다.
※질의서열을 직접 ftp로 다운받으려고 시도했지만 잘 안되는 관계로 아래 블로그에서 다운로드 받았다
http://bioprofiler.tistory.com/
ls -al #설치된 faa파일목록을 확인함.
ls *.faa > list_faa # 여러개의 파일목록을 lsit_faa에 담는다.
cat list_faa #목록을 확인한다.
, 이제, blastp를 사용하는 스크립트를 작성한다. 바로 blastp 명령문을 반복 실행하는 스크립트를 작성할 수도 있으나, 처음이고 하니 확인도 할 겸 아래의 스크립트를 mkblp.sh라는 이름의 파일로 (다른 이름도 상관 없다) ~/BLAST/query 디렉터리에 저장한다
1. vi mkblp_sh # mkblp_sh파일을 만든다. 내용은 아래처럼 입력한다. 입력할 때는 i, 나올때는 esc --> :wp를 입력
#!/usr/bin/env bash
genes=$(cat $1)
for gene in $genes
do
echo "blastp -query $gene -db ../blastdb/drosoph.aa -out ${gene%.faa}_blp.out"
done
2. 이제, 출력 결과를 확인하기 위해 실행해 보자. 이번에는 지금까지와 다르게 chmod 명령을 사용하지 않고 스크립트를 바로 실행하는 방법으로 해 본다.
sh mkblp.sh list_faa # sh명령어로 실행한다.
지금까지와 다르게 chmod 명령을 사용하지 않고 스크립트를 바로 실행하는 방법으로 해 본다.
list_faa는 mkblp.sh 스크립트의 $1에 해당하는 입력 매개 변수이다. 화면에 출력되는 내용을 보면, 다음과 같이 여러 줄의 blastp 명령문임을 알 수 있다.
$gene과 ${gene%.faa}_blp.out의 (변수) 위치에 list_faa 파일의 질의 서열 파일 이름이 들어가 있는 것을 알 수 있다. 이렇게 화면에 출력되는 내용을 파일로 저장하려면 >를 이용한 출력 방향 재지정 (redirection)을 사용하면 된다.
3. sh mkblp.sh list_faa > exeblp.sh # exeblp.sh는 잠시 전 화면에 출력되었던 내용을 저장한 파일이다.
cat 명령으로 내용에 이상이 없는 지 확인할 수 있다. 이렇게 작성된 텍스트 파일에는 앞 부분에 파일 처리에 필요한 처리기 지정 부분이 없다.
(편법일 수도 있지만) sh 명령으로 이 파일 (exeblp.sh) 역시 다음과 같이 실행할 수 있다.
이번에는 핵산 서열로 단백질 데이터베이스를 검색하여 그 결과를 html파일로 출력하는 예를 몇 가지 옵션을 지정해서 실행해 보자. 최근에는 mRNA (messenger RNA) 서열 데이터 (cDNA; complementary DNA)를 분석해야 하는 실험 방법들도 많은데, 이 경우에 활용할 수 있는 방법이기도 하다.
1. vi blxexe.sh # 텍스트파일을 작성한다.
중간에 스크립트를 파일로 만들지 않고 바로 명령을 실행하는 방법을 알아 보기로 하자. 다음의 내용을 ~/BLAST/query 디렉터리에 blxexe.sh이라는 이름의 텍스트 파일 (다른 이름도 상관 없다)로 저장한다.
2. cat blxexe.sh #준비된 핵산 질의 서열 파일을 확인하고 이들의 이름을 텍스트 파일로 저장하여 그 내용을 확인한다
3. sh blxexe.sh list_fna
blxexe.sh를 실행하여 for 루프로 blastx 명령문들을 실행한다. 잠시 기다려서 실행이 끝나면, 결과 html 파일을 확인한다.
여기서 바로 앞의 방식과 달라진 것은 for 루프 내에서 echo 명령으로 실행문을 출력하지 않고 바로 blastx 명령문을 실행한 것이다.
<에러발생>
sh blxexe.sh list_fna를 실행하는 과정에서 에러를 뿜는다. 왜그럴까? 구글링을 통해 이유를 찾는라고 1시간넘게 서치중이다.ㅠㅠ
. blxexe.sh 파일오타! #line6에 syntax오류라고 뜬다.
. #!/usr/bin/env bash
genes=$(cat $1)
for gene in $genes
do <-----------------요것을 빼먹은 엄청난 실수를 ㅠㅠ
blastx -query $gene -db ../blastdb/drosoph.aa -out ${gene%.fna}_e09_blx.html \ -evalue 1e-9 -html # do와 done 사이의 이 라인은 한 줄이다.
done ------
evalue와 le-9에 문제가 있는거 같다는 추측을 하게 된다.
결국 다음처럼 코드를 수정했다. -evalue le-9을 제거했다.
5. AK97891 단백질 서열을 nr 데이터베이스를 대상으로 E-value 옵션 1e-9로 검색한 html 결과의 일부이다.
C:\cygwin\home\user\BLAST\query로 찾아가서 html파일 하나를 더블클릭하면 다음처럼 문서가 나타난다.
'생명과학을 위한 생물학+정보학' 카테고리의 다른 글
| BLAST 데이터베이스_drosoph.aa 파일 (0) | 2018.03.09 |
|---|---|
| NCBI BLAST+ 데이터베이스만들기 (0) | 2018.03.09 |
| Cygwin_패키지설치 (0) | 2018.03.08 |
