FASTAおよびFASTQ形式のファイルの処理を、主にシェルのコマンドやsed,awk等のワンライナーを用いて行う。

1 ファイル名変換

fastqのファイル名は一様ではないが大抵の場合、サンプルIDなどの情報が_のような区切り文字で区切られている。ファイル名の文字列置換やcut,sedを使った切り出しでファイル名を作成する。
- array=(${変数/置換前文字列/ })でファイル名を分割した配列を作成する。
- sed 's///g'で文字列置換
- cutでファイル名を分割
sedは環境によって挙動が違うのでperlを使う方がトラブルが少ないかもしれない。

r1=C1_S244_L001_R1_001.fastq.gz

# Read1と対応するRead2のファイル名
echo ${r1/_R1/_R2}

# idを取り出す
echo $r1 | cut -f1,2 -d"_"
echo $r1 | perl -F"_" -lane 'print $F[0]'

# 拡張子削除 sed/perl
echo $r1 | sed -e 's/\..*//'
echo $r1 | perl -pe 's/^(.*?)\..*$/$1/g' # 最短マッチ
echo $r1 | perl -pe 's/^(.*)\..*$/$1/g'  # 最長マッチ(default)

# 拡張子削除 bashの変数内置換
echo ${r1%.*} 
echo ${r1%%.*}
echo ${r1#*.}
echo ${r1##*.}

2 複数のファイルをforループで処理

ペアエンドのfastqファイルを処理する場合、Read1のファイル名の配列を作って連続処理する。Read2はループの中で文字列置換してファイル名を作成して呼び出す。
${変数名/置換前文字列/置換後文字列},

# Read1のファイル名の配列を連続処理する。 
R1s=(*_R1_001.fastq.gz)
for r1 in ${R1s[@]} ; do
  r2=${r1/_R1/_R2}
  printf "%s\t%s\n" ${r1} ${r2}
done

# ファイル名配列の添字でループする場合
R1s=(*_R1_001.fastq.gz) ; R2s=(*_R2_001.fastq.gz)
for i in ${!R1s[@]} ; do
  r1=${R1s[$i]}
  r2=${R2s[$i]}
  echo -e "${r1}\\t${r2}"
done

# # ファイル一覧をループする場合の例
# for r1 in *_R1_001.fastq.gz; do echo $r1; done
# for r1 in `ls *_R1_001.fastq.gz`; do echo $r1; done
# for r1 in $(ls *_R1_001.fastq.gz); do echo $r1; done
# for r1 in $(ls ./ | grep "_R1_"); do echo $r1; done

## C1_S244_L001_R1_001.fastq.gz C1_S244_L001_R2_001.fastq.gz
## P1_S248_L001_R1_001.fastq.gz P1_S248_L001_R2_001.fastq.gz
## U1_S252_L001_R1_001.fastq.gz U1_S252_L001_R2_001.fastq.gz
## C1_S244_L001_R1_001.fastq.gz C1_S244_L001_R2_001.fastq.gz
## P1_S248_L001_R1_001.fastq.gz P1_S248_L001_R2_001.fastq.gz
## U1_S252_L001_R1_001.fastq.gz U1_S252_L001_R2_001.fastq.gz

# ファイル名とリード数
for r1 in *_R1_001.fastq.gz; do 
  r2=${r1/_R1/_R2}
  id=$(echo $r1 | cut -f1 -d"_")
  rn1=$((`gunzip -c $r1 | wc -l`/4))
  rn2=$((`gunzip -c $r2 | wc -l`/4))
  echo -e "${id}\\t${rn1}\\t${rn2}"
done

## C1   32961   32961
## P1   22480   22480
## U1   43996   43996

3 fastqデータ処理

3.1 fastqファイル

ls *.fastq.gz

## C1_S244_L001_R1_001.fastq.gz
## C1_S244_L001_R2_001.fastq.gz
## P1_S248_L001_R1_001.fastq.gz
## P1_S248_L001_R2_001.fastq.gz
## U1_S252_L001_R1_001.fastq.gz
## U1_S252_L001_R2_001.fastq.gz

3.2 fastqデータ処理

fastqファイルがgz圧縮されている場合、解凍結果をファイルではなく標準出力へ出力し、パイプで別コマンドに渡して処理 gzcat(Mac), zcat(linux)
awk '(NR-1)%4==1'で配列行だけを処理する。
fastqファイルの行数を4で割る計算をbashで行なう。
- echo "計算式" | bcで数値計算を行う。 echo $((計算式)), echo $[計算式]
- $()は、中の文字列をコマンドとして解釈し、その実行結果を返す。バッククォートで囲んでも良い。
- fastqがgzファイルではない場合wc -l file.fastqでは行数だけを渡せないのでcat file.fastq | wc -lとする。
awkでやる場合ENDブロックで計算するよりwcのほうが速い。
wc -l fileだと空行を数えてしまう。 grep . file | wc -l の方が良いかも。
bashでファイルを1行づつ処理する場合 while read line; do [処理]; done

# input
R1=C1_S244_L001_R1_001.fastq.gz

gunzip -c ${R1} | head -4
## gzcat ${fq} | head -4

# idを取得
gunzip -c ${R1} | awk '(NR-1)%4==0{print}' | head -3
## gunzip -c ${R1} | awk '(NR-1)%4==0' | head -3

# リードナンバー以降を除去
gunzip -c ${R1} | awk '(NR-1)%4==0{print $1}' | head -3

# リード長
gunzip -c ${R1} | awk '(NR-1)%4==1{print length($1)}' | head -3

# fastaに変換
gunzip -c ${R1} | awk '(NR - 1) % 4 < 2' | sed 's/^@/>/' | head -2

# リード数
echo $((`gunzip -c ${R1} | wc -l `/4))
echo $(gunzip -c ${R1} | wc -l)/4 | bc
## echo `gunzip -c ${R1} | wc -l`/4|bc
## echo $((`gunzip -c ${R1} | wc -l`/4))
## echo $(($(gunzip -c ${R1} | wc -l)/4))
## echo $[`gunzip -c ${R1} | wc -l`/4]
## gunzip -c ${R1} | wc -l | awk '{print $1/4}' 
## gunzip -c ${R1} | awk '{n++}END{print n/4}'


# 全てのリードについてリード長を取得する。
gunzip -c ${R1} | awk 'NR % 4 == 2 {print length($0)}'| head -n 3
## gunzip -c ${R1} | awk 'NR % 4 == 2' | head |  while read line; do echo ${#line}; done

3.3 fastqの先頭からサブセットを取り出す。

sed -n で対象行以外は出力しない(デフォルトでは処理しなかった行はそのまま出力)。pコマンドで処理した内容を出力する。
- sed -n '1,12p' 先頭から12行目まで(始めの3リード)を出力。
- sed 1,12d 例えば、1～12行目を削除する

for r1 in *_R1_001.fastq.gz
do
 r2=${r1/_R1/_R2}
 id=$(echo ${r1} | cut -f1 -d"_")
 r1sub=${id}_sub_R1.fastq
 r2sub=${id}_sub_R2.fastq
 gunzip -c ${r1} |sed -n '1,12p' > ${r1sub}
 gunzip -c ${r2} |sed -n '1,12p' > ${r2sub}
done

3.4 fastqからランダムにリードを取出す。

1. pasteでペアエンドのリードを横方向に結合
2. awkで4行ごとにタブで連結する。
3. shufでランダム行を抽出する。
4. リード1とリード2を別々に出力する。

<(gunzip -c fastq.gz)のようにプロセス置換を使って、解凍されたfastqファイルを作らないようにする。
awkのprintfは値をフォーマットして出力(改行がつかない)。行番号で分岐して1Reedの情報を tab区切りで出力する。
ランダム行を抽出はsort -Rでも良いと思う。sort -R | head -n 100
シード値固定でshufコマンド実行する場合--random-source=<(yes 123)
awkにシェルの変数を渡す場合
- awk -v 変数1=${変数1} -v 変数2=${変数2}
- x=123; echo "value:" | awk '{print $0 "'$x'"}'
awkの内側でパイプを介してファイル出力する場合{print | "gzip > outfile"}

# single endの場合
R1='C1_S244_L001_R1_001.fastq.gz'
R1SUB='C1_R1_sub.fastq.gz'
n=100
gunzip -c $R1 \
| awk '{ printf("%s",$0); n++; if(n%4==0) { printf("\n");} else { printf("\t");} }' \
| shuf -n $n \
| awk -F"\t" '{print $1"\n"$2"\n"$3"\n"$4 }' \
| gzip > $R1SUB

# paired-endの場合
R1='C1_S244_L001_R1_001.fastq.gz'
R2='C1_S244_L001_R2_001.fastq.gz'
R1SUB=R1_sub.fastq.gz
R2SUB=R2_sub.fastq.gz
n=100
RAND=123
paste <(gunzip -c $R1) <(gunzip -c $R2) \
| awk '{ printf("%s",$0); n++; if(n%4==0) { printf("\n");} else { printf("\t");} }' \
| shuf -n $n --random-source=<(yes ${RAND}) \
| awk -F"\t" '{print $1"\n"$3"\n"$5"\n"$7 | "gzip > '$R1SUB'"; \
print $2"\n"$4"\n"$6"\n"$8 | "gzip > '$R2SUB'"}'

3.5 fastqファイルからidを指定してリードを取り出す。

3.5.1 idのファイルを作成

# 取り出すリードのidファイル(id.txt)を作成。
# リード1,2いずれからも取り出す場合を考慮して、リードの情報以降を除去する。
R1='C1_S244_L001_R1_001.fastq.gz'
gunzip -c $R1 | awk '(NR-1)%4==0{print $1}' | sort -R | head -2 > id.txt
cat id.txt

## @M03099:22:000000000-APTMA:1:1101:14570:15532
## @M03099:22:000000000-APTMA:1:2109:13478:11843

3.5.2 fastqから指定idのリードを取得

awkで2つのファイルを比較する方法を利用する NR==FNR{1つめのファイル処理}NR!=FNR{2つめのファイル処理}

fastqファイルをタブ区切りに変換
idのファイルを読み込む際にidをkeyにしてハッシュを作成
tabに変換したfastqを読み込みながらfastqのid行でb.で作成したハッシュを検索
- リードの情報はスペースで区切られていることを想定している。
- tabに変換したfastqをawkで読み込む際にはtab区切りで読み込まないことに注意

@M03099:23:000000000-APTLR:1:1101:19252:1849 1:N:0:244
--> @M03099:23:000000000-APTLR:1:1101:19252:1849   この部分をidとして検索する

# fastqから指定idのリード抽出
R1='C1_S244_L001_R1_001.fastq.gz'
R2='C1_S244_L001_R2_001.fastq.gz'
ID='id.txt'

gunzip -c $R1 \
| awk '{ printf("%s",$0); n++; if(n%4==0) { printf("\n");} else { printf("\t");} }' \
| awk 'NR==FNR{a[$1]=$1}NR!=FNR{if ($1 in a) print $1" "$2"\n"$3"\n"$4"\n"$5; }' ${ID} -

gunzip -c $R2 \
| awk '{ printf("%s",$0); n++; if(n%4==0) { printf("\n");} else { printf("\t");} }' \
| awk 'NR==FNR{a[$1]=$1}NR!=FNR{if ($1 in a) print $1" "$2"\n"$3"\n"$4"\n"$5; }' ${ID} -

## @M03099:22:000000000-APTMA:1:1101:14570:15532 1:N:0:244
## CCTACGGGAGGCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATGCCGCGTGTGTGAAGAAGGTCTTCGGATTGTAAAGCACTTTAAGTTGGGAGGAAGGGCATTAACCTAATACGTTGGTGTTTTGACGTTACCGACAGAATAAGCACCGGCTAACTCTGTGCCAGCTGCCGCGGTAATACAGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTTTGTTAAGTTGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCAAA
## +
## 68BCCGGGCGGGGGGFGFGGFGF@EFGFGGFGGGGDGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGFGGGCFGGGGGCFFGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGDGGGGGGF@EGGFGFGCFGDGGGFFGCFGGGGGGGGGGGFGGGGFGGFGFFFFFGAFF9
## @M03099:22:000000000-APTMA:1:2109:13478:11843 1:N:0:244
## ACGGGCGGTGTGTACACATCTCCGAGCCCACGAGACGACCCCACGAGACGAGCCCACTAGACGGGTTCACGAGACGGGCCCACGAGCCGAGCCCACGCTACGAGCTCACGAGACGCGCCTACACAACGCCCCCACGAGCCGACCCTACGTGACGAGCCCACCCCATCCACCCACGTGACGGGGCCACCGGACTTTCCCATGGGAAAAATTCGCGTCATCAACCCCAGCCGACACGCTACTACATTAGCATTGCATAATACGATATGAATGCCCGTGCTTTTATCTTGTTTTCTCGTCCAC
## +
## CCC@<:FCGGG7FGFGA<E<CFF@7B7@F9@@C++CF7++8@FC+87@:@C:7CFFE=,9,:,,++,:E<,BFG7:4+@+48?:7@+++46>FF<>F+++@3+F7+><<,3588+*5:3+8,,8,18*5*41=<1*=*1***4*<22<*<<<<<1**18*/*/*****+2;*;CE8@5C****/*/8*12***00++3++*,*0+*2**3+1**2/**+++***22*/C****2***2*+++3+0+<+00+2<F+0+<:+**2*100+*+<<0*:2;)*+++++320++**+**(*2)19
## @M03099:22:000000000-APTMA:1:1101:14570:15532 2:N:0:244
## GACTACAGGGGTATCTAATCCTGTTTGCTCCCCACGCTTTCGCACCTCAGTGTCAGTATTAGTCCAGGTGGTCGCCTTCGCCACTGGTGTTCCTTCCTATATCTACGCATTTCACCGCTACACAGGAAATTCCACCACCCTCTACCATACTCTAGTCAGTCAGTTTTGAATGCAGTTCCCAGGTTGAGCCCGGGGATTTCACATCCAACTTAACAAACCACCTACGCGCGCTTTACGACCAGTAATTCCGATGAGCGCTTGCACCCTCTGTATTAGCGCTGCAGCGGTCACAGAGTTAGC
## +
## CBCCBFFFGGGG,CFFGDFGGGGFFGGGGGGGGDGGGGGGGGGGGGGGGFGGGGGGGGGGGGFG9EFFGGGGGGGGGEFFGGGEGFGFFFCEFCEGGGGGGGGGGGGCFCCCAFGGGGGGGGGGGGFGGGGGGGGFFGGGFGGGGGGGGGCGGFGGGGGFFGGGGGFD@FFGAFFGGCCGGGFGCF8FDFFGGGGCEE95=EEAEGG9EF,/>BFF>>7@FD>:A*5=9B;16:.>;(.2B(7:6.)7+;0;.*)08(4,3:23(4310:ABA+))(-((-65(.(((((46.4)003.)
## @M03099:22:000000000-APTMA:1:2109:13478:11843 2:N:0:244
## ACGCTATCTAACCATCTTCGATCCCCTGCCTCCCGTTCCTTCTTACTGCAGACATCCTTTGCACCTACCCTCTCGGTATTTACACTTCCGTAACTCCACGCATCTCGGCTCTTACACTCTACTCCTTGTACCCCACACGCCTCCTCTTAACTATTACTTCTCACCTCGCACTCCACATACTACAACCTCCCCTCCTATTTTATTCTCTCCTACTCCTATCTTCGACTACAGACATTCCTCCCGCCCTCTGTTTGTTTCCTACTAGCACTCGTATTTCCGTGCCTCGCCTCGCGTCTGTCTT
## +
## 6-8,6,,;CFE,;,;;E,@@C,,,,,C,,6;,,,+6,;6CEEA,CAFC,;,,;,C,6CC,,C,,,;,,,,,:,,,,8A,,:,,,:F9,6,9,,,:<@,,,48@+;,,,?D,?,,?CE,9+9:?,B?,,:,9,,A++4A+++,+04,,@84+6,,6+4,,,,,2,,6*+++++,,*3+*,@9A+,,83*+,@+=@:,+5,7,*,*;5,**)**3*2)**,32,5)<:)-.2)1))**4/4,6/>45@(())/*0)/,())))).))(/)0((/,(())+))/(())((,(,-(-)(**+)))

3.5.3 idのファイルとfastqファイルを引数として、指定idのfastqを取り出す関数

# idのファイルとfastqファイルを引数として、指定idのfastqを取り出す関数
function fqGet (){
    id=$1; fq=$2 
    cat ${fq} \
    | awk '{ printf("%s",$0); n++; if(n%4==0) { printf("\n");} else { printf("\t");} }' \
    | awk 'NR==FNR{a[$1]=$1}NR!=FNR{if ($1 in a) print $1" "$2"\n"$3"\n"$4"\n"$5; }' ${id} -
}

# 実際に実行する場合
R1='C1_S244_L001_R1_001.fastq.gz'
ID='id.txt'
gunzip -c ${R1} | fqGet ${ID} -

## @M03099:22:000000000-APTMA:1:1101:14570:15532 1:N:0:244
## CCTACGGGAGGCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATGCCGCGTGTGTGAAGAAGGTCTTCGGATTGTAAAGCACTTTAAGTTGGGAGGAAGGGCATTAACCTAATACGTTGGTGTTTTGACGTTACCGACAGAATAAGCACCGGCTAACTCTGTGCCAGCTGCCGCGGTAATACAGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTTTGTTAAGTTGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCAAA
## +
## 68BCCGGGCGGGGGGFGFGGFGF@EFGFGGFGGGGDGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGFGGGCFGGGGGCFFGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGDGGGGGGF@EGGFGFGCFGDGGGFFGCFGGGGGGGGGGGFGGGGFGGFGFFFFFGAFF9
## @M03099:22:000000000-APTMA:1:2109:13478:11843 1:N:0:244
## ACGGGCGGTGTGTACACATCTCCGAGCCCACGAGACGACCCCACGAGACGAGCCCACTAGACGGGTTCACGAGACGGGCCCACGAGCCGAGCCCACGCTACGAGCTCACGAGACGCGCCTACACAACGCCCCCACGAGCCGACCCTACGTGACGAGCCCACCCCATCCACCCACGTGACGGGGCCACCGGACTTTCCCATGGGAAAAATTCGCGTCATCAACCCCAGCCGACACGCTACTACATTAGCATTGCATAATACGATATGAATGCCCGTGCTTTTATCTTGTTTTCTCGTCCAC
## +
## CCC@<:FCGGG7FGFGA<E<CFF@7B7@F9@@C++CF7++8@FC+87@:@C:7CFFE=,9,:,,++,:E<,BFG7:4+@+48?:7@+++46>FF<>F+++@3+F7+><<,3588+*5:3+8,,8,18*5*41=<1*=*1***4*<22<*<<<<<1**18*/*/*****+2;*;CE8@5C****/*/8*12***00++3++*,*0+*2**3+1**2/**+++***22*/C****2***2*+++3+0+<+00+2<F+0+<:+**2*100+*+<<0*:2;)*+++++320++**+**(*2)19

3.5.4 idのアレイとfastqファイルを引数として、指定idのfastqを取り出す関数

# idをアレイとして持っている場合
function fqGet2 (){
    id=(${1}); fq=$2
    cat ${fq} \
    | awk '{ printf("%s",$0); n++; if(n%4==0) { printf("\n");} else { printf("\t");} }' \
    | awk 'NR==FNR{a[$1]=$1}NR!=FNR{if ($1 in a) print $1" "$2"\n"$3"\n"$4"\n"$5; }' \
    <(echo ${id[@]}|tr ' ' '\n') -
}
ids=($(cat ${ID}))
gunzip -c ${R1} | fqGet2 "${ids[*]}" -

3.6 不均一なペアエンドリードを同期させる

ペアエンドリードを個別にqcすると、ペア同士にならないリードができる。共通idのみからなるペアの配列を取り出す。
seqkitを使う場合は以下のようにする。
- seqkit common $r1 $r2 > $r1c
- seqkit common $r2 $r1c > $r2c

# idをアレイとfastqファイルを引数として、指定idのfastqを取り出す関数
function fqGet2 (){
    id=(${1}); fq=$2
    awk 'NR==FNR{a[$1]=$1}NR!=FNR{if ($1 in a) print $1" "$2"\n"$3"\n"$4"\n"$5; }' \
    <(echo ${id[@]}|tr ' ' '\n') \
    <(cat ${fq} | awk '{ printf("%s",$0); n++; if(n%4==0) { printf("\n");} else { printf("\t");} }')
}

# 不均一なペアエンドリードを同期させる
function readCom (){
    r1=$1;  r2=$2; pfx=$3

    # r1,r2 の共通id
    cid=($(cat <(cat ${r1} | awk '(NR-1)%4==0{print $1}') \
    <(cat ${r2} | awk '(NR-1)%4==0{print $1}') \
    | sort | uniq -c | awk '$1==2{print $2}'))

    # r1,r2のどちらかにしか存在しないidがある場合
    rid=($(cat <(cat ${r1} | awk '(NR-1)%4==0{print $1}') \
    <(cat ${r2} | awk '(NR-1)%4==0{print $1}') \
    | sort | uniq -c | awk '$1==1{print $2}'))

    if [[ ${#rid[@]} -ge 1 ]] ; then 
      fqGet2 "${cid[*]}" ${r1} > ${pfx}_c_R1.fq
      fqGet2 "${cid[*]}" ${r2} > ${pfx}_c_R2.fq
    else
    　echo -e '[ERROR] No need to pick common read.'
    fi
 
}

# input
r1='C1_S244_L001_R1_001.fastq.gz'
r2='C1_S244_L001_R2_001.fastq.gz'

# 不均一なペアエンドリードを作る
gunzip -c $r1 | sed -n '1,12p' > C1_qt_R1.fq 
gunzip -c $r2 | sed -n '5,16p' > C1_qt_R2.fq

# ペアのリードが揃っていないことを確認
echo "# Read1のID" ; awk '(NR-1)%4==0{print}' C1_qt_R1.fq
echo "# Read2のID" ; awk '(NR-1)%4==0{print}' C1_qt_R2.fq

readCom C1_qt_R1.fq C1_qt_R2.fq C1

# 共通IDのfastqが抽出されていることを確認
echo "# Read1のID(共通IDを抽出後)" ; awk '(NR-1)%4==0{print}' C1_c_R1.fq
echo "# Read2のID(共通IDを抽出後)" ; awk '(NR-1)%4==0{print}' C1_c_R2.fq

## # Read1のID
## @M03099:23:000000000-APTLR:1:1101:19252:1849 1:N:0:244
## @M03099:23:000000000-APTLR:1:1101:15851:1971 1:N:0:244
## @M03099:23:000000000-APTLR:1:1101:23102:1978 1:N:0:244
## # Read2のID
## @M03099:23:000000000-APTLR:1:1101:15851:1971 2:N:0:244
## @M03099:23:000000000-APTLR:1:1101:23102:1978 2:N:0:244
## @M03099:23:000000000-APTLR:1:1101:22637:2119 2:N:0:244
## # Read1のID(共通IDを抽出後)
## @M03099:23:000000000-APTLR:1:1101:15851:1971 1:N:0:244
## @M03099:23:000000000-APTLR:1:1101:23102:1978 1:N:0:244
## # Read2のID(共通IDを抽出後)
## @M03099:23:000000000-APTLR:1:1101:15851:1971 2:N:0:244
## @M03099:23:000000000-APTLR:1:1101:23102:1978 2:N:0:244

4 fasta形式のデータ処理

4.1 fastaファイル

4.2 fastaのヘッダ行を取り出して、文字列を分割する。

# fastaのヘッダ行を取り出して文字列分割
FA='C1_sub100.fasta'
grep "^>" ${FA} | tr -d '>' | head -3
grep "^>" ${FA} | tr -d '>' | cut -d ' ' -f1 | head -3
grep "^>" ${FA} | tr -d '>' | cut -d ' ' -f1 | cut -d ':' -f1,4,5,6,7 | head -3

## M03099:23:000000000-APTLR:1:1101:19252:1849 1:N:0:244
## M03099:23:000000000-APTLR:1:1101:15851:1971 1:N:0:244
## M03099:23:000000000-APTLR:1:1101:23102:1978 1:N:0:244
## M03099:23:000000000-APTLR:1:1101:19252:1849
## M03099:23:000000000-APTLR:1:1101:15851:1971
## M03099:23:000000000-APTLR:1:1101:23102:1978
## M03099:1:1101:19252:1849
## M03099:1:1101:15851:1971
## M03099:1:1101:23102:1978

4.3 fastaファイルの配列行に改行を追加する

fold半角文字列を指定した列数(バイト数)で改行させる。ヘッダー行の文字数が指定した数より多い場合、ヘッダー行も折りたたまれてしまう。
awkのsubstrを使って配列を指定数切り出すsubstr(対象文字列, 切り出し開始位置, 切り出す文字列の長さ
- awkでは -vオプションでコマンド実行前に変数を定義できる。

# foldを使う場合
FA='C1_sub100.fasta'
fold -w 60 ${FA} | head -3

# awkを使う場合
cat ${FA} | awk -v w=30 '/^>/{print $0}!/^>/{h=$0; while(h != 0){if(length(h) > w){print substr(h, 1, w); h=substr(h, w+1)}else{print h; h = 0} }}' | head -3

## >M03099:23:000000000-APTLR:1:1101:19252:1849 1:N:0:244
## CCTACGGGAGGCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATTC
## CCCGTGTGTGCTGCTTGCCTTCGGTTTGTATATCACTTTCTTTTTTGACGACTCTCTTTT
## >M03099:23:000000000-APTLR:1:1101:19252:1849 1:N:0:244
## CCTACGGGAGGCAGCAGTGGGGAATATTGG
## ACAATGGGCGAAAGCCTGATCCAGCCATTC

4.3.1 fastaファイルの配列部分の改行を削除

awk 'BEGIN{処理開始前の実行}条件式{処理1;処理2}END{最終行処理後の実行}'
awkのprintfは値をフォーマットして出力(改行がつかない)。 %d:整数値, %f:実数値, %s:文字列

# 元のfastaは配列部分の改行無しなので改行を追加する
FA='C1_sub100.fasta'
BRFA='C1_sub_br.fasta'
fold -60 ${FA} > ${BRFA}
head -3 ${BRFA}

# fastaの配列部分の改行を除く
cat ${BRFA} \
| awk '/^>/ { print n $0; n = "" }!/^>/ { printf "%s", $0; n = "\n" } END{ printf "%s", n }' \
| head -2

## >M03099:23:000000000-APTLR:1:1101:19252:1849 1:N:0:244
## CCTACGGGAGGCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATTC
## CCCGTGTGTGCTGCTTGCCTTCGGTTTGTATATCACTTTCTTTTTTGACGACTCTCTTTT
## >M03099:23:000000000-APTLR:1:1101:19252:1849 1:N:0:244
## CCTACGGGAGGCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATTCCCCGTGTGTGCTGCTTGCCTTCGGTTTGTATATCACTTTCTTTTTTGACGACTCTCTTTTTTTTATTACCTTTCTTTTTTTACTTTACCTACATTTTCTGCTCCGTCTTTCTTCTTTCCTTCTTCCTCTCTCATCCGTATGTTGCAACCTTTTTTCTTTTTTCCTTTTCTTTTATCTCTCTTTTTTGTTTCTTCACGTTTTCTTTTACCTCCCCTTTCTCTTCCTGTTAACTTCCTCCTAT

4.4 fastaファイルをtab区切りファイルに変換

# fastaファイル
FA='C1_sub100.fasta'

# fasta -> タブ区切りに変換
cat ${FA} | awk 'BEGIN{RS=">"; FS="\n"} NR>1 {print $1"\t"$2;}' | head -3

# 関数にする(fasta配列部分の改行削除 -> タブ区切りに変換 )
function faTab () {
  cat $1 \
  | awk '/^>/ { print n $0; n = "" }!/^>/ { printf "%s", $0; n = "\n" } END { printf "%s", n }' \
  | awk 'BEGIN{RS=">"; FS="\n"} NR>1 {print $1"\t"$2;}'  
}

## M03099:23:000000000-APTLR:1:1101:19252:1849 1:N:0:244    CCTACGGGAGGCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATTCCCCGTGTGTGCTGCTTGCCTTCGGTTTGTATATCACTTTCTTTTTTGACGACTCTCTTTTTTTTATTACCTTTCTTTTTTTACTTTACCTACATTTTCTGCTCCGTCTTTCTTCTTTCCTTCTTCCTCTCTCATCCGTATGTTGCAACCTTTTTTCTTTTTTCCTTTTCTTTTATCTCTCTTTTTTGTTTCTTCACGTTTTCTTTTACCTCCCCTTTCTCTTCCTGTTAACTTCCTCCTAT
## M03099:23:000000000-APTLR:1:1101:15851:1971 1:N:0:244    CCTACGGGAGGCAGCAGTGGGGAATATTGCACAATGGGCGCAAGCCTGATGCAGCCATGCCGCGTGTATGAAGAAGGCCTTCGGGTTGTAAAGTACTTTCAGCGGGGAGGAAGGTGTTTTGGTTAATAACCGCAGCAATTTACGTTACCCGCAGAAGATGCACCGGCTACCTCCGTGCCAGCAGCCGCGTTAATACGGAGGGTGCAAGCGTTAATCGGAATTACTGTGCGTAACGCGCTCGCAGTCGTTCTGTCAAGTCGGTTGTTAAATCCCCGGGCTCAACCTGGGAACTGCATTCTAA
## M03099:23:000000000-APTLR:1:1101:23102:1978 1:N:0:244    CCTACGGGAGGCAGCAGTGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGTGTGATGAAGGCCCTAGGGTTGTAAAGCTCTTTCCCCGGATAAGATAATGACGTTATCCGTAGAAGAAGCCCCTGCTAACTTCGTGCCAGCATCCGCGGTAATACGACTGGCGCTAGCGTTGTTCGGTATTACTGGTCGTAAAGCTCACGTAGGCTGATCGATCAGTCATGGGTGCAATCCCATGTCTCACCCCTGGCTCTGCCTTTCATACTTTCCCTCTGGATTATGTCAGAG

4.5 塩基配列操作

大文字小文字変換
相補鎖に変換

# revcomp 
function revcomp(){ echo $1 | tr ATGCYRKMWBVDHNatgcyrkmwbvdhn TACGRYMKWVBHDNtacgrymkwvbhdn | rev ; }
# faUpper
function faUpper (){ echo $1 | tr atgcyrkmwbvdhn ATGCYRKMWBVDHN ;}

4.6 GC/AT/N の割合

# fastaファイル
FA='C1_sub100.fasta'

# faGCAT, contents of [ GC | AT | N ]
function faNGCAT () {
  cat $1 \
  | awk '/^>/ { print n $0; n = "" }!/^>/ { printf "%s", $0; n = "\n" } END { printf "%s", n }' \
  | awk 'BEGIN{RS=">"; FS="\n"} NR>1 {print $1"\t"$2;}'  \
  | awk -F"\t" '{ALL=length($2); NAT=gsub(/[AT]/,"", $2); NGC=gsub(/[GC]/,"",$2); NN=gsub(/[Nn]/, "", $2);\
  NATGC=NAT+NGC; if(NATGC==0){GCC=0; ATC=0;}else{GCC=NGC/NATGC ; ATC=NAT/NATGC}; \
  print $1"\t"ALL"\t"ATC"\t"GCC"\t"NN/ALL }'
}
  
faNGCAT ${FA} | head -3

## M03099:23:000000000-APTLR:1:1101:19252:1849 1:N:0:244    301 0.598007    0.401993    0
## M03099:23:000000000-APTLR:1:1101:15851:1971 1:N:0:244    301 0.461794    0.538206    0
## M03099:23:000000000-APTLR:1:1101:23102:1978 1:N:0:244    301 0.468439    0.531561    0

4.7 multi-fastaファイルを指定配列数で分割する

sprintf(書式, 値1, 値2)
- sprintf("%s_%0" dgt "d" "." sfx, pfx, n )
- sprintf("%s_%03d.fasta, "C1_sub100", 1 ) # 変数展開して書くとこんな感じ. => C1_sub100_001.fasta

# faChunk
function faChunk () {
   size=$1; input=$2
   pfx=${input%%.*}
   sfx=${input##*.}
   dgt=$((`grep -c "^>" ${input} | wc -c`-1))
   cat ${input} \
   | awk -v size=${size} -v pfx=${pfx} -v dgt=${dgt} -v sfx=${sfx} \
   '/^>/ { n++; if (n % size == 1) { close(fname); fname=sprintf("%s_%0" dgt "d" "." sfx, pfx, n )}} \
   {print >> fname}'
}

# 20配列ずつ分割
FA='C1_sub100.fasta'
N=20
faChunk $N $FA
grep -c "^>" C1_sub100*

## C1_sub100.fasta:100
## C1_sub100_001.fasta:20
## C1_sub100_021.fasta:20
## C1_sub100_041.fasta:20
## C1_sub100_061.fasta:20
## C1_sub100_081.fasta:20

4.8 マルチファスタファイルからidで配列を検索して取り出す

# complete match of idlist and header of fasta 
function faGet (){
  id=$1; fa=$2 
  cat ${fa} \
  | awk '/^>/ { print n $0; n = "" }!/^>/ { printf "%s", $0; n = "\n" } END{ printf "%s", n }' \
  | awk 'BEGIN{RS=">"; FS="\n"} NR>1 {print $1"\t"$2;}' \
  | awk 'NR==FNR{a[$1]=$1}NR!=FNR{if ($1 in a) print ">"$1"\n"$2}' ${id} -
}

function faGet2 () {
  ids=(${1}); fa=$2
  cat ${fa} \
  | awk '/^>/ { print n $0; n = "" }!/^>/ { printf "%s", $0; n = "\n" }END{ printf "%s", n }' \
  | awk 'BEGIN{RS=">"; FS="\n"} NR>1 {print $1"\t"$2;}' \
  | awk 'NR==FNR{a[$1]=$1}NR!=FNR{if ($1 in a) print ">"$1"\n"$3}' <(echo ${ids[@]}|tr ' ' '\n') -
}


# idとマッチしたfastaとマッチしないfastaに分割 faGetrest <(cat ${ids[@]}) input.fasta
function faGetrest (){
  id=$1; fa=$2; hit=$3; rest=$4 
  cat ${fa} \
  | awk '/^>/ { print n $0; n = "" }!/^>/ { printf "%s", $0; n = "\n" } END{ printf "%s", n }' \
  | awk 'BEGIN{RS=">"; FS="\n"} NR>1 {print $1"\t"$2;}' \
  | awk 'NR==FNR{a[$1]=$1}NR!=FNR{if ($1 in a) print ">"$1"\n"$2 > "'${hit}'"; else print ">"$1"\n"$2 > "'${rest}'"}' ${id} -
}

fa='C1_sub100.fasta'
ids=(`cat ${fa} |grep ">"|sed -e "s/^>//" | cut -f1 -d" " | sed -n '3,5p'`)

function faGet2 () {
  ids=(${1}); fa=$2
  cat ${fa} \
  | awk '/^>/ { print n $0; n = "" }!/^>/ { printf "%s", $0; n = "\n" }END{ printf "%s", n }' \
  | awk 'BEGIN{RS=">"; FS="\n"} NR>1 {print $1"\t"$2;}' \
  | awk 'NR==FNR{a[$1]=$1}NR!=FNR{if ($1 in a) print ">"$1"\n"$3}' <(echo ${ids[@]}|tr ' ' '\n') -
}

faGet2 "${ids[*]}" ${fa}

## >M03099:23:000000000-APTLR:1:1101:23102:1978
## CCTACGGGAGGCAGCAGTGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGTGTGATGAAGGCCCTAGGGTTGTAAAGCTCTTTCCCCGGATAAGATAATGACGTTATCCGTAGAAGAAGCCCCTGCTAACTTCGTGCCAGCATCCGCGGTAATACGACTGGCGCTAGCGTTGTTCGGTATTACTGGTCGTAAAGCTCACGTAGGCTGATCGATCAGTCATGGGTGCAATCCCATGTCTCACCCCTGGCTCTGCCTTTCATACTTTCCCTCTGGATTATGTCAGAG
## >M03099:23:000000000-APTLR:1:1101:22637:2119
## CCTACGGGAGGCAGCAGTGGGGAATATTGCACAATGGGGGAAACCCTGATGCAGCCATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTCAGTTGTGAGGAAAAGTTAGTAGTTAATACCTGCTAGCCGTGACGTTAACAACAGAAGAAGCACCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGCTAGATGTGAAAGCCCTGGGCTCAACCTGGGATGGTCATTTAGA
## >M03099:23:000000000-APTLR:1:1101:18265:2244
## CCTACGGGAGGCAGCAGTGGGGAATATTGCACAATGGGCGCAAGCCTGATGCAGCCATGCCGCGTGTATGAAGAAGGCCTTCGGGTTGTAAAGTACTTTCAGCGGGGAGGAAGGTGTTGTGGTTAATAACCGCAGCAATTGACGTTACCCGCAGAAGAAGCACCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAA

FASTQ and FASTA file handling

S.Konishi

2023-02-09

1 ファイル名変換

2 複数のファイルをforループで処理

3 fastqデータ処理

3.1 fastqファイル

3.2 fastqデータ処理

3.3 fastqの先頭からサブセットを取り出す。

3.4 fastqからランダムにリードを取出す。

3.5 fastqファイルからidを指定してリードを取り出す。

3.5.1 idのファイルを作成

3.5.2 fastqから指定idのリードを取得

3.5.3 idのファイルとfastqファイルを引数として、指定idのfastqを取り出す関数

3.5.4 idのアレイとfastqファイルを引数として、指定idのfastqを取り出す関数

3.6 不均一なペアエンドリードを同期させる

4 fasta形式のデータ処理

4.1 fastaファイル

4.2 fastaのヘッダ行を取り出して、文字列を分割する。

4.3 fastaファイルの配列行に改行を追加する

4.3.1 fastaファイルの配列部分の改行を削除

4.4 fastaファイルをtab区切りファイルに変換

4.5 塩基配列操作

4.6 GC/AT/N の割合

4.7 multi-fastaファイルを指定配列数で分割する

4.8 マルチファスタファイルからidで配列を検索して取り出す