Shell script로 Fasta file 만들기

왜 해야됨..

https://elemental-tabletop-1a6.notion.site/2022_MutClust-b293d6b27c9545dfb762fe00c2bad643
학교에서 학부연구생으로 지금 진행하고 있는 프로젝트가 있는데, 질병청에서 받아온 데이터가 잘못되어 있는데가 있어서 데이터를 새로 만들기로 했다. 애초에 실험을 통해서 받은 fastq raw data가 있었기 때문에 실험이 이상하게 진행되지 않은 이상 이 데이터에는 문제가 없다.(없었으면 좋겠다..)

# fastq file format
      
1 @NB551490:285:HMJ3TAFX2:1:11101:3762:1095 1:N:0:CCATCATTAG+NGAGGCAACC
2 AAATTGATCTCCAGGCGGTGGTTTAGCACTAACTCTGGAAAAATCTGTATTATTAGGTGTATCAACATAACCTGTAGGTACAGCAACTAGGTTAACACCTGTAGAAAAACCTAGCTGTAAAGGTAAATTGGTACCAACAGCTTCTCTAG
3 +
4 AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEAEEEEEEEEEEEEE<EEEAAAAAEAA<EEEAAAAAAEAA
5 @NB551490:285:HMJ3TAFX2:1:11101:15573:1100 1:N:0:CCATCATTAG+AGAGGCAACC
6 CCCAAGAATAGCATAGATGCCTTCAAACTCAACATTAAATTGATGGGTGTTGGTGGCAAACCTTGTATCAAAGTAGCCACTGTACAGTCTAAAATGTCAGATGTAAAGTGAACATCAGTAGTCTTACTCTCAGTTTAGCAACAACTCAC
7 +
8 AAAAA/AAEEEEEAE6EAAE///EEEEEEEEE<EEEEEE6EA/6//<EA66/6AAA/EEE/AEA<EEAAEEE/EE/EAEEEEAE/E//6AAEEEAEE/E/AA/AAEAEA</EEEE/EEEE/A6/<AAA<////6/A/<6<A<6EAA6</

데이터에 대해서 간단하게 말을 하자면, 1번, 5번 줄의 정보는 NGS에서 얻은 header정보로 각 read에 대한 정보라고 보면 된다. 2번, 6번 줄은 염기서열에 대한 정보이고, 4번, 8번 줄은 각 염기서열의 정확도에 대한 정보이다. 이 정확도는 ASCII로 이루어져있어서 높을수록 좋다고 보면 된다. 이 데이터를 그래도 쓸 수가 없는 이유는 이 서열이 어디부분에 위치해 있는지 나타나있지 않기 때문에 alignment라는 과정을 통해서 위치를 찾고 그 위치에 맞게 재배치 해주는 과정을 해줘야 했다. 쉽게 말하면, 나는 이 데이터를 내가 변이를 쉽게 볼수 있도록 나타내기 위해 fasta 형태로 만들어야 했다.
아래에 보이는 NC_045512.2가 우리가 알고 있는 코로나 바이러스의 코드 번호이다.

# fasta file format
>NC_045512.2
NNTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAA
CGAACTTTAAAATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAAC
TAATTACTGTCGTTGACAGGACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTG
TTGCAGCCGATCATCAGCACATCTAGGTTTTGTCCGGGTGTGACCGAAAGGTAAGATGGAGAGCCTTGTC
CCTGGTTTCAACGAGAAAACACACGTCCAACTCAGTTTGCCTGTTTTACAGGTTCGCGACGTGCTCGTAC
GTGGCTTTGGAGACTCCGTGGAGGAGGTCTTATCAGAGGCACGTCAACATCTTAAAGATGGCACTTGTGG
CTTAGTAGAAGTTGAAAAAGGCGTTTTGCCTCAACTTGAACAGCCCTATGTGTTCATCAAACGTTCGGAT
...

Volume

애초에 이전에 했던 프로젝트가 fastq파일을 통해서 의미를 찾아내야하는 프로젝트였기 때문에, 위에서 해야하는 과정은 어렵지 않았다. trimmomatic이나 bwa와 같은 alignment를 위한 tool가 잘 만들어져 있기에 그 tool을 가져다 쓰기만 하면 됐기 때문이다. 근데 문제는 분석해야할 파일이 너무 많은것이었다. 환자 1명에 대한 데이터의 개수가 8개인데 환자의 수가 300명이었기 때문에 대략 2400개의 데이터를 다뤄야했을 뿐만 아니라 데이터 하나당 용량이 100MB정도라 데이터 하나하나를 보면서 관리하는 것을 거의 불가능에 가까웠다. 처음에 자신있는 파이썬으로 과정을 진행해보려고 했지만 데이터 양뿐만 아니라 하나 당 용량도 커서 파이썬으로 하기에는 비효율적이었다.
연구실 서버가 linux커널이었기 때문에 bash를 이용해서 한번도 써본적 없는 shell script로 자동화를 해볼려고 했다.

Data preprocessing

생물데이터는 특이하게 paired-end라고 2개의 데이터를 1쌍으로 취급해서 데이터를 다룬다. 쉽게 말하면 그냥 input에 넣는 데이터가 1쌍이라고 생각하면 된다. 그래서 밑에 있는거 처럼 2개의 데이터를 한번에 넣어줘야 하는 과정을 shell script로 나타내야했다.

	COV-CCO-0011_COVIDseq_S3_L001_R1_001.fastq
	COV-CCO-0011_COVIDseq_S3_L001_R2_001.fastq

일단 fastq 데이터를 만든 NGS가 사용하는 adapter라는 것이 있는데, adapter가 뭐냐면 데이터를 만들기 위해서 임의로 붙인 염기서열이다. 이거는 covid에서 나올 수 없는 서열이기 때문에 쓰레기 데이터라고 볼 수 있다. 그래서 이 쓰레기 데이터를 없애주기 위해서 Trimmomatic이라는 도구를 사용해서 데이터를 전처리 했다.

	java -jar $trimmo/trimmomatic-0.39.jar PE -phred33 $left $right 
    $trimpath/$file/$trimleft.trimmed.fastq 
    $trash/$trimleft.up.trimmed.fastq 
    $trimpath/$file/$trimright.trimmed.fastq 
    $trash/$trimright.up.trimmed.fastq 
    ILLUMINACLIP:$trimmo/adapters/TruSeq3-PE.fa:2:30:10:2:True LEADING:3 

위에서 보는것처럼 $left $right 변수 2개가 필요한데 shell에서 for문을 쓸때는 데이터 하나밖에 가지고 오는 것 밖에 안되기 때문에 코드를 수정해서 2개씩 들어가도록 해야했다.

for fn in $samp_path/*.fastq
do
  file=`basename $fn`
  mkdir $trimpath/$file
  cnt=0
  for samples in $fn/*.fastq
  do
    temp=$((cnt%2))
    if [ $temp -eq 0 ]
    then
      left=$samples
    else
      right=$samples
    fi
    if [ $temp -eq 1 ] && [ $cnt -gt 0 ]
    then
      trimleft=`basename $left .fastq`
      trimright=`basename $right .fastq`
      java -jar $trimmo/trimmomatic-0.39.jar PE -phred33 $left $right $trimpath/$file/$trimleft.trimmed.fastq $trash/$trimleft.up.trimmed.fastq $trimpath/$file/$trimright.trimmed.fastq $trash/$trimright.up.trimmed.fastq ILLUMINACLIP:$trimmo/adapters/TruSeq3-PE.fa:2:30:10:2:True LEADING:3 TRAILING:3 MINLEN:36
    fi
    cnt=$((cnt+1))
    echo $trim
	done
done

for문을 이용해서 .fastq 형식의 데이터를 모두 가져온 후, cnt가 짝수이면 왼쪽, cnt가 홀수이면 오른쪽으로 가게 한 이후에 이 과정을 한번 거칠때 마다 preprocessing을 해주는 과정을 통해서 자동화를 진행했다. 가장 오래 걸렸던 부분은 if안의 구절에 띄어쓰기를 제대로 해주지 않으면 코드가 진행되지 않는 것을 찾는데 오래 걸렸다.(bash할 때는 무조건 띄어쓰기부터 고려 해줘야 될 듯..)
참고로 basename은 파일 경로에서 파일이름만 가져오는 것인데, 나는 새로운 파일을 만들어줘야 하는 작업이기 때문에 이름을 그래도 옮기기 위해서 사용했다.

Alignment

Alignment 하는 과정도 파일 2개를 하나로 합치는 과정이기 때문에 input이 2개가 필요했다. 그래서 마찬가지로 위의 코드와 유사하게 진행을 했다. alignment를 하는 tool을 bwa인데 앞에서 말했듯이 이 서열의 위치가 어디에 있는지를 알려주는 tool이라고 보면 된다.

//sam information format
<QNAME> <FLAG> <RNAME> <POS> <MAPQ> <RNEXT> <CIGAR> <PNEXT> <TLEN> <SEQ> <QUAL> <TAG>
NB551490:285:HMJ3TAFX2:1:11106:18529:15900
99      NC_045512.2     709     60      149M    =       804     244     
GCTTGGCACTGATCCTTATGAAGATTTTCAAGAAAACTGGAACACTAAACATAGCAGTGGTGTTACCCGTGAACT
CATGCGTGAGCTTAACGGAGGGGCATACACTCGCTATGTCGATAACAACTTCTGTGGCCCTGATGGCTACCCTC
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEE
EEEEEEEEEEEEEEEEEEAEEEEEEEE/EEEEEEEEEAEE<EEEEEEEAAEA<<<EAEE<EEEEEEEAA<AEE/  
NM:i:0  MD:Z:149        MC:Z:149M       AS:i:149        XS:i:0

쉽게 설명하자면 709번 위치부터 149 길이의 서열이라는 것을 알려주는 것이라고 생각하면 된다. 다른 정보들도 알 수 있지만 나는 위치만 알면 되기 때문에 이 내용만 언급하도록 하겠다.

for fn in $trimpath/*.fastq
do
  file=`basename $fn .fastq`
  mkdir $mappingpath/$file
  mkdir /data3/projects/2020_MUTCLUST/FASTQ/fasta/$file
  cnt=0
  for samples in $fn/*.fastq
  do
    temp=$((cnt%2))
    if [ $temp -eq 0 ]
    then
      left=$samples
    else
      right=$samples
    fi
    if [ $temp -eq 1 ] && [ $cnt -gt 0 ]
    then
      trimleft=`basename $left .fastq`
      trimright=`basename $right .fastq`
      echo $trimleft
      echo $trimright
      $bwa mem -t 20 -v 2 $bwaidx $left $right | $samtools sort -o $mappingpath/$file/$trimleft.bwa.bam
    fi
    cnt=$((cnt+1))
    #echo $trim
	done
done

코드가 다른게 있다면 proprocessing한 데이터를 가져와야했기 때문에 거기에 맞는 dirname을 써준거 밖에 없다. 그리고 bwa라는 alignment tool을 이용해서 sam file을 생성하고

$bwa mem -t 20 -v 2 $bwaidx $left $right

sam은 용량이 너무크기 때문에 압축을 하여 데이터를 저장하는데, 그 데이터 형식이 bam이다. 그래서 samtools를 이용해서 sam을 bam으로 변환해주는 작업도 진행했다.

$samtools sort -o $mappingpath/$file/$trimleft.bwa.bam

fasta file 생성

결국 내가 만들고 싶은 데이터는 fasta의 형태이기 때문에 bam을 fasta로 만드는 작업을 진행해야했다. 다행히 samtools에 concensus라는 옵션으로 코드를 작성하면 fasta file을 만들어준다는 것을 알아서 그것을 사용해보기로 했다.

samtools consensus -a --show-ins no $files -o $path.fasta

$files는 bam이고 결과는 .fasta로 나오게 된다. 마지막으로 다뤄야 했던 데이터인지라 폴더명을 최대한 간단하게 해서 나타내고 싶어서 최대한 중복되는 파일명을 지운뒤에 작업을 진행했다.

for fn in $bamdir/*
do
  bampos=`basename $fn 1_COVIDseq`
  mkdir $fastadir/$bampos
  for files in $fn/*.bam
  do
    filename=`basename $files _R1_001.trimmed.bwa.bam`
    path=$fastadir$bampos/$filename
    echo $path
    samtools consensus -a --show-ins no $files -o $path.fasta
  done  
done

Result

>NC_045512.2
NNTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAA
CGAACTTTAAAATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAAC
TAATTACTGTCGTTGACAGGACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTG
TTGCAGCCGATCATCAGCACATCTAGGTTTTGTCCGGGTGTGACCGAAAGGTAAGATGGAGAGCCTTGTC
CCTGGTTTCAACGAGAAAACACACGTCCAACTCAGTTTGCCTGTTTTACAGGTTCGCGACGTGCTCGTAC
GTGGCTTTGGAGACTCCGTGGAGGAGGTCTTATCAGAGGCACGTCAACATCTTAAAGATGGCACTTGTGG
CTTAGTAGAAGTTGAAAAAGGCGTTTTGCCTCAACTTGAACAGCCCTATGTGTTCATCAAACGTTCGGAT
GCTCGAACTGCACCTCATGGTCATGTTATGGTTGAGCTGGTAGCAGAACTCGAAGGCATTCAGTACGGTC
GTAGTGGTGAGACACTTGGTGTCCTTGTCCCTCATGTGGGCGAAATACCAGTAGCTTACCGCAAGGTTCT
...

이런 식으로 데이터를 얻을 수 있었다. 다행히 데이터가 많았지만 쓰레기값이나 오류는 없어서 중간에 끊기지 않고 진행할 수 있었다. 덕분에 bash와 shell script로 작업하는 경험을 해볼 수 있어서 나중에 있을 데이터 관리에 대한 자신감을 기를수 있었을뿐만 아니라 파이썬에 대한 의존성도 덜 수 있었다.

추가 작업

위에 데이터에서 보는 것처럼 N으로 이루어진 mapping이 되지 않은 region이 있는데, 이것을 어떻게 해주어야할지 고민을 해봐야겠다.