google colabolatoryでcircRNAを同定したい④

# circtoolsは結局断念
circtoolsに関してはこちらのページを参照にしました。

https://pypi.org/project/circtools/

このページには、circtoolsについてこう書かれていました。
>circtoolsは、データ処理コードのほとんどをPython 3 (>=3.7)で、プロットと統計解析の部分はR (>=4.0.0)で書かれている。このツールは多くの外部依存性を持っており、そのほとんどが標準的なバイオインフォマティクスツールとパッケージである。インストールはデフォルトで、必要な依存関係をすべてインストールする。

一つ引っかかる。「ほとんどをpython3」:thinking:
一部は違うということか:thinking::thinking:
githubのサイトを見てみると、こう書いてある。

https://github.com/dieterich-lab/circtools?tab=readme-ov-file

>circtoolsパッケージはPython3(>=3.4)で書かれており、detectとreconstructの2つのモジュールはPython2(>=2.7)のインストールが必要である。

え、またか:weary:
うーん、やめとこう。
## CIRI2を使う

初心に戻って、これまで使われてきているpipelineでpython3ベースのものはないか探すことにしました。（なぜ最初からそうしなかったのか）
ということで、気を取り直して、**CIRI2**を使用することにしました。

https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0571-3#Sec18

上記の原著論文によると
>CIRIはPerlで実装され、LinuxとMac OS Xで広範囲にテストされている。CIRIを使うのに他のツールは必要ない。CIRIはCIRI-simulatorとともにパッケージ化されており、[34]で入手可能である。

とのことで、これは「エラーを避けられるのではないか？」と期待が高まります:flushed:

このpipelineの原著論文のタイトルは上記のように 
**CIRI: an efficient and unbiased algorithm for *de novo* circular RNA identification** 
と一見、de novoでcircRNAを同定できそうな雰囲気がありますが、**Materials and Methods**の**Algorithm**では、

>CIRIは、FASTAフォーマットされた参照配列と、BWA-MEMアルゴリズム[17]によって生成されたSAMアライメントの2種類のファイルを必要とする。

と書かれています。おそらくは新規に同定するという意味で*de novo*と表記しているようで、 ***de novo* assembly** 自体は自分でやらないといけないようです。

**CIRI2**は`CIRI-cookbook`として懇切丁寧な取扱説明サイトがありました。(もちろん、他のpipelineでも形は違えど、チュートリアルがあると思います。)

https://ciri-cookbook.readthedocs.io/en/latest/CIRI2.html#an-example-of-running-ciri2

こちらのサイトの**An example of running CIRI2**を真似てみることにしました。

## CIRI2のインストール
まずCIRI2をgithubを介してinstallしました。
```bash:クローン化
!git clone https://github.com/bioinfo-biols/CIRI-full.git
```
クローン化したフォルダは`/content/drive/MyDrive`の直下に置きました。
## sampleデータのダウンロード
次に`MyDrive`の階層にsampleフォルダを作りました。
(ただ、わざわざここで作らなくても、`サンプルデータのダウンロード`のコード内で`os.makedirs(download_dir, exist_ok=True)`を加えて、ファイルの保存先ディレクトリが存在しない場合作成するとしてもよかったですね。)
```python:sampleフォルダの作成
import os
print(os.getcwd())
os.chdir('/content/drive/MyDrive/sample')
```
sampleデータはCIRI2のフォルダ内には含まれておらず、別のサイトからダウンロード`https://sourceforge.net/projects/ciri/files/CIRI2/`する必要がありました。

>Please download CIRI2 and test_data2.zip from https://sourceforge.net/projects/ciri/files/CIRI2/ and then unzip the three files in test_data2.zip.

test_data2.zipをダウンロードしていきます。`.sraファイル`と異なり`wget`ではダウンロードが出来なかったので、chatGPT先生に教えてもらって、先ほど作ったsampleフォルダにダウンロードしました。初心者の私には`status_code`が新鮮です。
```python:サンプルデータのダウンロード
import requests
#import os 
#先にosをimportしていたので、今回は割愛

# ダウンロードするファイルのURL
url = "https://sourceforge.net/projects/ciri/files/CIRI2/CIRI_v2.0.6.zip/download"

# ダウンロードするファイルの保存先ディレクトリ
download_dir = "/content/drive/MyDrive/sample"

# ダウンロードするファイルの保存先パス
file_path = os.path.join(download_dir, "CIRI_v2.0.6.zip")

# ファイルをダウンロード
response = requests.get(url)

# レスポンスが成功した場合
if response.status_code == 200:
    # ダウンロードしたファイルを保存
    with open(file_path, 'wb') as f:
        f.write(response.content)
    print("ファイルのダウンロードが完了しました。保存先:", file_path)
else:
    print("ファイルのダウンロードに失敗しました。")
```
sampleフォルダ内にダウンロードできているか確認します。
```bash:sampleフォルダ内の確認
!ls 
#カレントディレクトリ以外のディレクトリの中身を知りたい場合は
#!ls $path
```
`CIRI_v2.0.6.zip`が表記されていたので、無事ダウンロードできていることが確認できました。
次に`unzip`で解凍します。

```bash:zipの解凍
!unzip test_data2.zip 
#引数として -d $pathでoutput folderを指定できる
!ls
```
`chr1.fa`   `CIRI_v2.0.6.zip`  `test.sam`
`chr1.gtf`  `CIRIerror.log`     `__MACOSX`
と表示されました
各ファイルの内容としては、こんな感じでしょうか。

`chr1.fa：1番染色体のfastaファイル`
`test.sam：サンプルの配列をマッピングしたデータ`
`chr1.gtf：chr1をアノテーションしたデータ`
`__MACOSX：zip解凍時に作られるフォルダ。windowsユーザーには関係ないみたい。`
`CIRIerror.log：エラーのlog`

次にCIRI2を使用していきます。
## CIRI2で解析
CIRI2の使い方は以下の通り。
```bash:How to run CIRI2
!perl CIRI2.pl -I in.sam -O outfile -F ref.fa (-R ref_dir/)
```
引数は以下の通り。
>-I, --in
    input SAM file name (required; generated by BWA-MEM)
-O, --out
    output circRNA list name (required)
-F, --ref_file
    FASTA file of all reference sequences. Please make sure this file is 
    the same one provided to BWA-MEM. Either this argument or -R/--ref-dir 
    is required.
-R, --ref_dir
    directory of reference sequence(s). Please make sure fasta files in 
    this directory are from the FASTA file(s) provided to BWA-MEM. Either 
    this argument or -F/--ref-file is required.
-A, --anno
    input GTF/GFF3 formatted annotation file name (optional)
    
CIRI2自体はperlで書かれているためか、path通しがうまくできず、だらだらと絶対パスを使用して起動させました。(アノテーションファイルの指定が無くても解析可能なようですが、今回はアノテーションファイルを指定する方法で行いました。)

```bash:circRNAの同定
!perl /content/drive/MyDrive/CIRI-full/bin/CIRI_v2.0.6/CIRI2.pl -I test.sam -O outfile -F chr1.fa -A chr1.gtf
```
`[Tue Apr  9 03:42:48 2024] CIRI finished its work. Please see output file outfile for detail.`
の表記が出たので再度sampleフォルダの中身を確認しました。
`outfile`と`outfile.log`が増えていました。この`outfile`の中身は以下の情報が入っています。ちなみに`tab`で区切られています。

<table>
  <caption>outfile</caption>
  <thead>
    <tr>
      <th>circRNA_ID</th><th>chr</th><th>circRNA_start</th><th>circRNA_end</th><th>#junction_reads</th><th>SM_MS_SMS</th><th>#non_junction_reads</th><th>junction_reads_ratio</th><th>circRNA_type</th><th>gene_id</th><th>strand</th><th>junction_reads_ID</th>
    </tr>
  </thead>
  <tr>
    <td> chr1:16770127|16775696 </td> <td>chr1</td><td>16770127</td><td>16775696</td><td>8</td><td>5_6_0</td><td>192</td><td>0.077</td><td>exon</td><td>ENSG00000157191.15</td><td>+</td><td>GFGG-GA2-1:2:33:2000:1585,GFGG-GA2-1:2:33:2000:1587...</td>
  </tr>
  <tr>
  <td> chr1:16775588|16778510 </td> <td>chr1</td><td>16775588</td><td>16778510</td><td>6</td><td>5_5_0</td><td>192</td><td>0.056</td><td>exon</td><td>ENSG00000157191.15</td><td>+</td><td>GFGG-GA2-1:2:33:2000:1276,GFGG-GA2-1:2:33:2000:1277...</td>
  </tr>
</table>

`outfile`には拡張子がありませんが、おそらく`.txtファイル`か`.csvファイル`形式のように思います。
どちらのファイル形式として保存しても、開くことができました。

ちなみに各項目については<a href ="https://ciri-cookbook.readthedocs.io/en/latest/CIRI2.html#an-example-of-running-ciri2">cookbook</a>に以下のような説明が載っていました。
<table><tr><th>列</th><th>説明</th></tr>
<tr><td>1</td><td>ID of a predicted circRNA in the pattern of "chr:start</td></tr>
<tr><td>2</td><td>chromosome of a predicted circRNA</td></tr>
<tr><td>3</td><td>start loci of a predicted circRNA on the chromosome</td></tr>
<tr><td>4</td><td>end loci of a predicted circRNA on the chromosome</td></tr>
<tr><td>5</td><td>circular junction read (also called as back-spliced junction read) count of a predicted circRNA</td></tr>
<tr><td>6</td><td>unique CIGAR types of a predicted circRNA. For example, a circRNAs have three junction reads: read A (80M20S, 80S20M), read B (80M20S, 80S20M), read C (40M60S, 40S30M30S, 70S30M), then its has two SM types (80S20M, 70S30M), two MS types (80M20S, 70M30S) and one SMS type (40S30M30S). Thus its SM_MS_SMS should be 2_2_1.</td></tr>
<tr><td>7</td><td>non-junction read count of a predicted circRNA that mapped across the circular junction but consistent with linear RNA instead of being back-spliced</td></tr>
<tr><td>8</td><td>ratio of circular junction reads calculated by 2#junction_reads/(2#junction_reads+#non_junction_reads). #junction_reads is multiplied by two because a junction read is generated from two ends of circular junction but only counted once while a non-junction read is from one end. It has to be mentioned that the non-junction reads are still possibly from another larger circRNA, so the junction_reads_ratio based on it may be an inaccurate estimation of relative expression of the circRNA.</td></tr>
<tr><td>9</td><td>type of a circRNA according to positions of its two ends on chromosome (exon, intron or intergenic_region; only available when annotation file is provided)</td></tr>
<tr><td>10</td><td>ID of the gene(s) where an exonic or intronic circRNA locates</td></tr>
<tr><td>11</td><td>strand info of a predicted circRNAs (new in CIRI2)</td></tr>
<tr><td>12</td><td>all of the circular junction read IDs (split by ",")</td></tr></table>

おー、うまくいった。テストデータなので当然かもしれませんが、今まではテストですらうまくいきませんでした。大きく前進した気持ちになります。

拡張子を書き足したときに`.txtファイル`として保存してみたので、このデータを`.csvファイル`へ書き出したいと思います。
```python:.csvファイルへの書き出し
import pandas as pd

# データを読み込む
data = pd.read_csv('outfile.txt', sep='\t')

# .CSVファイルに書き出す
data.to_csv('output.csv', index=False)
```
カレントディレクトリである`sampleフォルダ`内に無事`outfile.csv`として書き出されました。