google colabolatoryでcircRNAを同定したい②（step1~4）

前回は、circRNAFinderのgitをクローンするところで終わりました。

https://qiita.com/sakiimai01/items/d650741f1fd23a9c60e5

今回はいよいよcircRNAFinder内のTutrial.shの指示に従い、一つずつステップを踏んでいきたいと思います。基本的には**材料集めとそのプロセシング**の回になります。
# Tutorial.shを眺める
Tutrialとしてどんなことが書かれているのか見ていきたいと思います。
```bash
#! /usr/bin/bash
# This tutorial helps you get started with circRNAFinder by demonstrating how to detect, annotate and visualize circRNAs using public RNA-Seq data. It also could be used as a template to build a pipeline for your own data.

# Step1: Create a working directory
･･･
# Step2: Download and build human genome/transcript sequences and gene annotations. (You only have to do it once.)
･･･
## Step3: Create fasta folder, download the public RNA-Seq data from SRA database.
･･･
## Step4: Convert .sra files into FASTA format, collapse the identical reads.
･･･
## Step5: Edit configure files
･･･
## Step6: Run circRNAFinder
･･･
```
げげっ・・・:fearful:
step 6まである・・・:rolling_eyes:
しかも、circRNAFinderを起動させるまでの処理が長い:sob:
でも、頑張るしかあるまい:nerd:
## Step1:ワーキングディレクトリを作る
ここはなんとかできそう。
まずカレントディレクトリを確認します。
```python:python
#カレントディレクトリを確認
import os
print(os.getcwd())
```
`/content`にあることがわかりました。
続けて、元々提供されていたコードに従い、`example`というフォルダを`/content/drive/ MyDrive`の中に作りたいと思います。ですので、まず、MyDriveをワーキングディレクトリとしたいと思います。
```python:python
# MyDriveディレクトリへ移動
os.chdir('/content/drive/MyDrive')
```
```bash:bash
#新しくディレクトリを作成
!mkdir example
!cd example
```
```python:python
#新しく作ったディレクトリに移動
os.chdir('/content/drive/MyDrive/example')
```
## Step2:ヒトゲノム/トランスクリプトの配列及び遺伝子アノテーションをダウンロードしビルドする。
ここでも指示通り、先ほど作った`example`フォルダの下層に`build`フォルダを作成しました。
```bash:bash
!mkdir build
```
再度、ワーキングディレクトリを移動します。
```python:python
os.chdir ('/content/drive/MyDrive/example/build')
```
次にpath通しをして、git内にある`SeqAnnoDownloadBuild.py`を使用します。
```python:python
#path通し
import sys
sys.path.append('/content/drive/MyDrive/circRNAFinder/src')
```
```bash:bash
# genome/transcript、gene annotationをダウンロード&ビルド
!SeqAnnoDownloadBuild.py hg19 --download --genome
!SeqAnnoDownloadBuild.py hg19 --build --genome
!SeqAnnoDownloadBuild.py hg19 --download --transcriptome
!SeqAnnoDownloadBuild.py hg19 --build --transcriptome
!SeqAnnoDownloadBuild.py hg19 --download --annotation
!SeqAnnoDownloadBuild.py hg19 --build --annotation
!cd ..
```
ここで、何度かエラーがでました。
どうも`SeqAnnoDownloadBuild.py`内に書かれている`print`の指示の文法が古いことが原因のようで、`print`から`print()`に変更し保存した後、実行するとうまくいきました。同様のエラーが他の`.pyファイル`でも起きています。
`build`フォルダ内に`hg19.gtf`、`hg19.gtf.build`、`hg19.grf.gz`、`hg19.tar.gz`、`hg19_dna.fa`、`hg19_trans.fa`というファイルができていました。
## Step3: fastaフォルダを作り、そこにSRAデータベースで公開されているRNA-Seq dataをダウンロードする
先ほど同様に`build`フォルダ内に`fasta`フォルダを作成し、移動します。
念のため、移動した後、ディレクトリの位置を確認します。
```python:python
print(getcwd())
```
`/content/drive/MyDrive/example/build/fasta`となっており、移動が確認できました。
指示には、「マニュアルでデータをダウンロードするか、もしAsperaをインストールしているなら、次のコマンドでデータをダウンロードできる」と書いてありましたが、コードが複雑であることと、普通に`wget`でデータをダウンロードできることから、指示には従わず次のコードを実行しました。(自動でプログラムやファイル名が表示されたのでTabキーを押しまくり、ゴリゴリなコードとなりました。)
```bash:bash
# Or if you've installed the Aspera, you can download them using the following commands in terminal:
# ~/.aspera/connect/bin/ascp -QTr -k1 -l640M -L log --retry-timeout=10 -i ~/.aspera/connect/etc/asperaweb_id_dsa.putty anonftp@ftp-private.ncbi.nlm.nih.gov:sra/sra-instant/reads/ByRun/sra/SRR/SRR650/SRR650317/SRR650317.sra .

#ここからが私のコード。圧縮した形式でダウンロード。
!fastq-dump --split-files SRR650317 --gzip
!fastq-dump --split-files SRR364679 --gzip
!fastq-dump --split-files SRR364680 --gzip
!fastq-dump --split-files SRR364681 --gzip
!fastq-dump --split-files SRR384963 --gzip
!fastq-dump --split-files SRR384964 --gzip
!fastq-dump --split-files SRR384965 --gzip
```
`fasta`内に`SRR650317`～`SRR384965`の`fastq.gz`がダウンロードされていました。
(尚、pair-endのものは`SRR650317`のみで、これは`SRR650317_1.fastq.gz`と`SRR650317_2.fastaq.gz`に分かれ、他は、single-endのため、例えば`SRR364679_1.fastq.gz`のように一つだけファイルができていました。)
## Step4: Convert .sra files into FASTA format, collapse the identical reads
Step4では`sra2fasta.py`と`fastx_collapser.py`を使用するよう指示がありましたが、結局使いませんでした。
`sra2fasta.py`は「`sra.file → fasta.file`」へ変換するプログラムですが、`fastq-dump`がその機能を担ってくれているので、使う必要はありませんでした。
`fastx_collapser.py`は、複数の圧縮された`fasta.gz`または`fastq.gz`ファイルを取得し、解凍して一つの`fasta.file`にまとめるプログラムでした。
しかし、python2ベースのpipelineだからか、エラーが頻発し、先ほどの`print`の記述を変えたり、`args = parser.parse_args()`を`args = parser.parse_args(args=[])`に変換したりしてみましたが、エラーがでてうまくいきませんでした。
```python:エラーの内容
usage: colab_kernel_launcher.py [-h] [-f | -q] input output
colab_kernel_launcher.py: error: the following arguments are required: input, output
An exception has occurred, use %tb to see the full traceback.

SystemExit: 2
```
そこで解凍、`seqtk`を用いての`fastaファイル`への変換、マージは一つ一つ段階を追って行うことにしました。`seqtk`は`.gzファイル`にも対応しているので、解凍する必要はなかったな、と後から気づきました。
```bash:bash
#ディレクトリ内のファイルを全て解凍
!gzip -d *.fastq.gz

#seqtkでfastq→fastaへ変換

#seqtkをinstall
!conda install -c bioconda seqtk -y

#フォルダ内の.fastqを.fastaへ変換(以降chatGPT作)
import os

# フォルダのパス
folder_path = '/content/drive/MyDrive/example/build/fasta'

# フォルダ内のファイルを取得
files = os.listdir(folder_path)

# フォルダ内の各ファイルについて処理
for file in files:
    # .fastqファイルのみ処理
    if file.endswith('.fastq'):
        # 入力ファイルと出力ファイルのパスを作成
        input_file = os.path.join(folder_path, file)
        output_file = os.path.join(folder_path, os.path.splitext(file)[0] + '.fa')
        
        # seqtkを使って変換
        cmd = f'seqtk seq -a {input_file} > {output_file}'
        print(f'Converting {input_file} to {output_file}')
        os.system(cmd)
```

```bash:bash
#各fastqを一つにまとめる
# HEK293 (Paired-end reads) 
%cat SRR650317_1.fa SRR650317_2.fa > HEK293.fa

# CD19 (Single-end reads) #
%cat SRR364679.fa SRR384963.fa > CD19.fa

# CD34 (Single-end reads) #
%cat SRR364680.fa SRR384964.fa > CD34.fa

# neutrophils (Single-end reads) #
%cat SRR364681.fa SRR384965.fa > neutrophils.fa
%cd ..
```
ここまでで一通り、材料集めは終わりました。次回、いよいよcircRNAFinderを起動しようと思います。