cutadapt_demo

---

---

cutadaptインストール

---

• 以下のサイトを参考にしてインストール(https://cutadapt.readthedocs.io/en/stable/installation.html) macでは$HOME/Library/Python/3.9/bin にインストールされた。

# pipを使ってインストール
python3 -m pip install --user --upgrade cutadapt

# バージョン指定してインストール
pip install cutadapt==1.18

# 確認
cutadapt --version
cutadapt --help

アダプターの位置を確認

---

fastqファイル

• PRJEB19901のデータ(アンプリコンシーケンスリード, Miseq, 16SrRNA, V3-V4)を使用

リードの端の配列を確認

• リードの端にprimer配列があるかをgrepで確認。
  • 341F-CCTACGGGNGGCNGCAG
  • 515F-GTGYCAGCMGCCGCGGTAA
  • 806R-GACTACNNGGGTATCTAATCC
• 混合塩基の部分は正規表現で記述　[ATGC]\{2\}
• プライマーに5’位置にランダム配列(0~5塩基)が挿入される場合[ATGC]\{0,5\}
• grepの結果をパイプでmoreとかに渡すと色が消えてしまうので--color=alwaysをつける。
• リードスルーを考慮する場合複数箇所をgrepする　grep -e "xx" -e "yy"

# input
R1='SRR5446818_1.fastq.gz'
R2='SRR5446818_2.fastq.gz'

# read1、read2の5'末端にプライマーがあるか確認する
gunzip -c ${R1} | grep -B1 "CCTACGGG[ATGC]GGC[ATGC]GCAG" --color=always | head
gunzip -c ${R2} | grep -B1 "GACTAC[ATGC]\{2\}GGGTATCTAAT" --color=always | head

# forward primerにランダム配列(0~5塩基)が含まれる
gunzip -c ${R1} | grep -B1 "[ATGC]\{0,5\}CCTACGGG[ATGC]GGC[ATGC]GCAG" --color=always | head

# read1の5末にreverseプライマー/read2の5末にforwardプライマーが含まれることがある。インデックスホッピング?
gunzip -c ${R1} | grep -B 1 "GACTAC[ATGC]\{2\}GGGTATCTAATCC" --color=always | head
gunzip -c ${R2} | grep -B 1 "CCTACGGG[ATGC]GGC[ATGC]GCAG" --color=always | head

# リードスルーを考慮してforward/reverseプライマーの２箇所を検索する
gunzip -c ${R1} | grep -B1 -e "GGATTAGATACCC[ATGC]\{2\}GTAGTC" \
| grep -B1 -e "CCTACGGGAGGCAGCAG" -e "GGATTAGATACCC[ATGC]\{2\}GTAGTC" \
--color=always | head

gunzip -c $R2 | grep -B1 -e "CTGCTGCCTCCCGTAGG" \
| grep -B1 -e "GACTAC[ATGC]\{2\}GGGTATCTAAT" -e "CTGCTGCCTCCCGTAGG" \
--color=always

cutadapt実行

• アダプタ配列を指定。5末(-g)と3末(-a)でcutの仕方が違う。両方のパターンでcutする場合は-bを使う。
  • -g, -G 5’末端のアダプターの配列(-g:リード1, -G: リード2)。
  • -a, -A 3’末端のアダプター配列(-a:リード1, -A:リード2)。アダプター以降の塩基はトリミングされる。
  • -b, -B 5’または3’末端のアダプターの配列。-aと-gの両方のタイプのマッチが可能(-b:リード1, -B: リード2)
  • -o, -p トリミングされたリード(-o:リード1, -p: リード2)
    
• -n: 各リードから除くアダプターの個数[default:1]
  
• ペアエンドで実行する場合-a ${FORWARD}...${RCREVERSE} -A ${REVERSE}...${RCFORWARD}
  • この場合--discard-untrimmedでプライマーを含まない配列を削除
• 古いバージョンでは上記のようなペアエンドアンプリコンのプライマー除去が使えない(<=1.18?)。その場合-g, -aを両方を指定する。ただしデフォルトではよりマッチする方の1つが除かれる為、-n 2を指定する必要がある。cutadapt -n 2 -g ${FORWARD} -a ${REVERSE} ${R1} > ${R1CUT}
• ただし、場合によっては3末のアダプターが２回除かれることが起こりうる。その場合パイプで２回cutadaptを実行する。cutadatp -g ${FORWARD} $R1 | cutadapt -a ${RCREVERSE} - > ${R1CUT}
  • https://cutadapt.readthedocs.io/en/stable/recipes.html
• アンカー付きアダプターについてヘルプを参照すること

# 相補鎖配列をプリントする関数
function revcomp () { 
  echo $1 | tr ATGCYRKMWBVDHNatgcyrkmwbvdhn TACGRYMKWVBHDNtacgrymkwvbhdn | rev 
}

# Primer配列
FP='CCTACGGGNGGCNGCAG'
#FP='NNNNNCCTACGGGNGGCNGCAG'
RP='GACTACNNGGGTATCTAAT'
FP3p=$(revcomp ${FP})
RP3p=$(revcomp ${RP})

# input/output
R1='SRR5446818_1.fastq.gz'
R2='SRR5446818_2.fastq.gz'
R1CUT=${R1/.fastq.gz/_cut.fastq.gz}
R2CUT=${R2/.fastq.gz/_cut.fastq.gz}

# paired-endで実行
cutadapt -j 4 -a ${FP}...${RP3p} -A ${RP}...${FP3p} --discard-untrimmed -o ${R1CUT} -p ${R2CUT} ${R1} ${R2} > log1.txt

# single endで実行(これだと例えば２回reverseプライマーがカットされることが起こりうる。)
cutadapt -n 2 -g ${FP} -a ${RP3p} -o ${R1CUT} ${R1}
cutadapt -n 2 -g ${RP} -a ${FP3p} -o ${R2CUT} ${R2}

# single endで5末、3末のカットを連続で実行(--discard-untrimmedをつけるとペアの不均一ができる)
cutadapt -g ${FP} ${R1} | cutadapt -Z -a ${RP3p} -o ${R1CUT} - >>log3.txt
cutadapt -g ${RP} ${R2} | cutadapt -Z -a ${FP3p} -o ${R2CUT} - >> log3.txt


# single endで5末、3末のカットを連続で実行(--discard-untrimmedをつけるとペアの不均一ができる)
cutadapt -j 4 -g ${FP} ${R1} 2>log3.txt| cutadapt -j 4 -Z -a ${RP3p} -o ${R1CUT} - >>log3.txt
cutadapt -j 4 -g ${RP} ${R2} 2>>log3.txt| cutadapt -j 4 -Z -a ${FP3p} -o ${R2CUT} - >> log3.txt


# passしたリード数を確認
echo $((`gunzip -c $R1|wc -l`/4)); echo $((`gunzip -c $R1CUT|wc -l`/4))

FP='CCTACGGGNGGCNGCAG'; RP='GACTACNNGGGTATCTAAT'
FP3p=$(revcomp ${FP}); RP3p=$(revcomp ${RP})

for R1 in ./fastq/*_1.fastq.gz; do 
R2=${R1/_1.fastq.gz/_2.fastq.gz}
R1CUT=`basename ${R1}| sed -e "s/.fastq.gz/_cut.fastq.gz/"`; 
R2CUT=`basename ${R2}| sed -e "s/.fastq.gz/_cut.fastq.gz/"`;
ID=`basename ${R1} | sed -e "s/_1.fastq.gz//"`
echo -e $R1"¥t"$R2"¥t"${ID}"¥t"$R1CUT"¥t"$R2CUT
cutadapt -j 4 --discard-untrimmed \
-a ${FP}...${RP3p} -A ${RP}...${FP3p} \
-o ${R1CUT} -p ${R2CUT} \
${R1} ${R2} > ${ID}_cut.log
done