[Cutadapt] 02. Cutadapt로 adapter를 자르고 결과를 간단하게 정리하기

수정 : 2023.05.08

내용 추가: 2024.01.10 

 

▶ 이전 글 : Cutadapt설치하기

[Cutadapt] 01.Linux에 Cutadapt 4.2 ver 설치하기

16s rRNA 데이터를 분석하기 위해서 시퀀싱 과정에서 서열에 부착한 adapter서열을 잘라내야 합니다. 이를 위해 가장 많이 사용되는 프로그램이 Cutadapt입니다. 최신 버전은 4.2로 Linux와 Window, macOS 도 지원합니다. 그러나 대부분 많은 분석 tool이 Linux환경에서 구동되기 때문에 command line과 반드시 친해져야만 합니다. 

 

이제 cutadapt를 사용한 시퀀싱 파일의 adapter를 제거하고 제거된 파일의 총 몇 %가 남았는지 나타내는 결과 파일도 정리해 보겠습니다. 

 

 

🟦 Cutadapt를 쓰기 전 

illumina사의 시퀀서를 예로 우리에게 필요한 것은 아래와 같습니다.

 

1. demultiplexing 된 fastq파일

- multiplexing은 시퀀서에 같이 돌리는 서열들을 구분하기 위해 붙여둔 분자 바코드를 붙이는 단계입니다. 이러한 분자 바코드를 이용해 서열을 나누어 주는 단계를 Demultiplexing이라 부릅니다. 보통 시퀀싱 회사에 서열을 맡기면 Demultiplexing이 끝난 sample.fastq.gz파일을 전송해 줍니다. 

 

2. fastq파일과 원래 샘플의 이름 리스트

- name.txt파일로 저장해 리눅스로 옮겨 놉니다.

 

3. adapter 서열

- 이는 각 kit에 따라, 원하는 영역에 따라 달라집니다.

- 이 게시글의 예로 bacteria 16s rRNA의 V1-V3 region영역을 시퀀싱 할 때에는 

      8F : 5'-AGAGTTTGATCCTGGCTCAG-3'

      534R : 5'-ATTACCGCGGCTGCTGG-3’

      의 primer서열을 사용하였습니다.

 

이제 adapter 서열을 기존 fastq서열에서 잘라내 주도록 하겠습니다.

 

🟦 Cutadapt로 adapter 자르기 

사용할 파일을 아래와 같다

 

1.  fastq.gz 형태로 되어있는 시퀀스 

001AcCUSw_1_1.fastq.gz
001AcCUSw_1_2.fastq.gz
001AcLUSw_3_1.fastq.gz
001AcLUSw_3_2.fastq.gz
...

2. fastq.gz 파일들의  이름이 적힌 txt파일(namelist.txt)

#tr -d '\r' < name.txt > name.txt

$ less name.txt

001AcLUSw_3
001AcRUSw_1
..

 

 

일단 분석이 사용되는 곳에 위치를 만들어 줍니다.

mkdir fastq/
mkdir Cutadapt/fastq

cd name.txt Cutadapt/

sample의 이름이 적힌 파일을 Cutadapt폴더로 이동시킵니다.

우리가 adapter를 자르고자 하는 파일은 Cutadapt/fastq/안에 이동시킵니다.

 

 

이제 Cutadapt를 사용해 봅시다! 일단 illumina 사의 시퀀서를 이용할 때는 한 샘플당 두 개의 서열이 나오기 때문에 총 4개의 primer를 잘라내야 합니다. 

# Cutadapt 실행
source cutadapt-venv/bin/activate

for i in  *_1.fastq.gz; do  \
  SAMPLE=$(echo ${i} | sed "s/_1.fastq.gz//g") ;   \
  echo "On sample : $SAMPLE" ; \  # 어떤 샘플의 adapter를 자르고있는지 실시간으로 표시하기 위해
  cutadapt  \
  -g ^AGAGTTTGATCCTGGCTCAG \ 
  -G ^ATTACCGCGGCTGCTGG   \
  -a CCAGCAGCCGCGGTAAT \
  -A CTGAGCCAGGATCAAACTCT \
  --minimum-length 20   \
  -Q 20 -q 20  \  # forward와 reverse 서열의 QC 기준을 20으로 잡음
  -o  Cutadapt/fastq/${SAMPLE}_1.fastq.gz  \
  --paired-output  Cutadapt/fastq/${SAMPLE}_2.fastq.gz   \
  ${SAMPLE}_1.fastq.gz  ${SAMPLE}_2.fastq.gz   \
  >>  Cutadapt/result.txt 2>&1 ;  done  # 결과 요약파일 저장

각 값을 살펴보면 아래와 같다.

 

  -g ^AGAGTTTGATCCTGGCTCAG  # 잘라낼 forward read의 5' adapter
  -G ^ATTACCGCGGCTGCTGG   # 잘라낼 reverse read의 5'adapter

  -a CCAGCAGCCGCGGTAAT # 잘라낼 forward read의 3'adapter
  -A CTGAGCCAGGATCAAACTCT # 잘라낼 reverse read의 3'adapter
  --minimum-length 20   # read길이 20이 안되면 버린다 (나중에 dada2 분석 시 error를 막기 위함)
  -Q 20 -q 20   # forward와 reverse 서열의 QC 기준을 20으로 잡음
  -o ~   # forward read 저장 위치 
  --paired-output  ~ # reverse read 저장 위치

 

어뎁터에 대한 자세한 옵션값은 링크를 참고하길 바랍니다. (https://cutadapt.readthedocs.io/en/stable/guide.html)

 

 

이제 result.txt파일을 봅시다. 

less Cutadapt/result.txt

 

 

 

🟦 결과 한눈에 보기

Cutadapt에는 pass와 filtering단계가 있습니다. 

위의 result.txt파일을 이용해 최종 파일이 몇 %가 남았는지 확인해 봅시다.

일단 위의 result파일에서 빨간색 박스 안에 있는 값을 가져오려고 합니다. 이는 아래와 같은 스크립트를 사용하여 각 값을 가져오고 히우 합쳐줍니다. 이때 파일의 name.txt가 있어야 파일을 깔끔하게 정리할 수 있습니다.

 

paste name.txt <(grep "passing filters" cutadapt_01.txt | cut -f3 -d "(" | tr -d ")") <(grep "filtered" cutadapt_01.txt | cut -f3 -d "(" | tr -d ")") > Output.txt
less -S Output.txt

 

less -S Output.txt

짜잔! 2열의 데이터를 보면 결과적으로 Cutadapt 이후 남은 샘플의 양을 확인할 수 있습니다.  

 

 

 

 

🟦 자세한 결과 요약하기 (ChatGPT)

# 현재 샘플 이름을 추출하여 변수에 저장
sample_names=$(ls *_1.fastq.gz | cut -f1 -d "_")

# 결과 파일 생성
{ echo -e "SampleID\ttotal_readpairs\tpass_read\tpass_per\ttotal_basepairs\tfilt_base\tfilt_per"; \
  paste <(printf "${sample_names}\n") \
        <(grep "Total read pairs processed" ./result.txt | awk '{print $5}') \
        <(grep "passing filters" ./result.txt | awk '{print $5}') \
        <(grep "passing filters" ./result.txt | cut -f3 -d "(" | tr -d ")") \
        <(grep "Total basepairs processed" ./result.txt | awk '{print $4}') \
        <(grep "filtered" ./result.txt | awk '{print $4}') \
        <(grep "filtered" ./result.txt | cut -f3 -d "(" | tr -d ")"); } > output.txt

 

 

위를 함수 형태로 만들어 보면 아래와 같습니다. 

#!/bin/bash

# 사용법 출력 함수
usage() {
    echo "Usage: $0 -s sample_names -p path -i input -o output"
    echo "  -s: Sample names file (e.g., names.txt)"
    echo "  -p: Path to the directory containing result.txt"
    echo "  -i: Input file (result.txt)"
    echo "  -o: Output file name"
    exit 1
}

# 옵션 파싱
while getopts "s:p:i:o:" opt; do
  case $opt in
    s) sample_names="$OPTARG" ;;
    p) path="$OPTARG" ;;
    i) input="$OPTARG" ;;
    o) output="$OPTARG" ;;
    \?) usage ;;
  esac
done

# 필수 옵션 확인
if [ -z "$sample_names" ] || [ -z "$path" ] || [ -z "$input" ] || [ -z "$output" ]; then
    usage
fi

# 샘플 이름을 변수에 저장
sample_names_content=$(cat "$sample_names")

# 결과 파일 생성
{
  echo -e "SampleID\ttotal_readpairs\tpass_read\tpass_per\ttotal_basepairs\tfilt_base\tfilt_per";
  paste <(printf "%s\n" "${sample_names_content}") \
        <(grep "Total read pairs processed" "${path}/${input}" | awk '{print $5}') \
        <(grep "passing filters" "${path}/${input}" | awk '{print $5}') \
        <(grep "passing filters" "${path}/${input}" | cut -f3 -d "(" | tr -d ")") \
        <(grep "Total basepairs processed" "${path}/${input}" | awk '{print $4}') \
        <(grep "filtered" "${path}/${input}" | awk '{print $4}') \
        <(grep "filtered" "${path}/${input}" | cut -f3 -d "(" | tr -d ")");
} > "${output}"

echo "Result file created: ${output}"

 

 

실행은 bash Cutadapt_result.sh -s sample_names -p ./ -i result.txt -o output2.txt 로 수행할 수 있습니다. 

 

---

🟦 Reference

- cutadapt 공식 가이드 : https://cutadapt.readthedocs.io/en/stable/guide.html

- https://www.biostars.org/p/254099/

- https://www.biostars.org/p/414888/

- https://astrobiomike.github.io/amplicon/dada2_workflow_ex