Python演習

目的

　バイオインフォマティクスのための大量情報処理においてはPythonのようなスクリプト言語を利用して、複数のプログラムで連続的に処理することがよくなされる。また、生命情報の多くの部分を占める配列情報の効率的な取り扱いにおいても、スクリプト言語は非常に強力なプログラミング言語である。
　そこで本実験では、以下の項目毎に1問ずつ課題を与え、実際にスクリプト言語によるプログラムの作成を行う。課題は、基礎的な6課題（必修）と応用的な3課題（希望者のみ）に分かれる。それぞれの課題を解く事によって、単純にスクリプト言語の習得だけでなく、バイオインフォマティクス研究を行うための基礎知識として、配列ファイルのフォーマットやいくつかの配列データベースの使い方、配列の類似性検索プログラムの利用法、タンパク質立体構造データベースからのデータの取得、統計量の計算など、バイオインフォマティクス研究に必要な基本的な知識の習得も目的とする。
（本課題は東京大学医科学研究所機能解析イン・シリコ分野で利用されている教育用課題を発展させたものである）

基礎的課題

▼　１ 基礎統計量の計算（スクリプト言語の基礎の基礎） 
標準入力から複数の整数を入力として受け取り、入力された整数の合計・平均値・標準偏差を計算し、標準出力へ書き出すプログラムを作成せよ。
入力は1行に1つの整数が与えられる。
プログラムの出力（合計、平均、標準偏差）は下のフォーマットに従うこと。

-- ここから
sum [タブ] [合計値] [改行]
mean [タブ] [平均値] [改行]
sd [タブ] [標準偏差] [改行]
-- ここまで

合計、平均、標準偏差以外の出力をしたいとき（入力値など）は標準エラー出力に書き出してください。
合計、平均、標準偏差以外を標準出力へ書き出すと、自動チェックスクリプトが正しく動きません。
また、なんらかの理由で値が計算できないときは [合計値] [平均値]  [標準偏差] のかわりに "NA"と表示してください。
実行例1
  • 入力: 1, 2, 3, 4, 5
  • 合計 = 15
  • 平均値 = 3
  • 標準偏差 = 1.581
実行例2
  • 入力: -2, -1, +1, +2
  • 合計 = 0
  • 平均値 = 0
  • 標準偏差 = 1.826
入力として不適切な例
  • 1.5
  • 1.0
  • 1+2
  • XVIII
▼　２ 配列操作の基礎（相補配列） 
マルチFASTA フォーマットの塩基配列データファイルのパスを引数にとり、その相補配列を5'→3'方向に出力するプログラムを作れ。
出力はマルチFASTAフォーマットとしてパースできる形にし、標準出力へ書き出すこと。
出力に含まれる配列ヘッダ（">"から始まる行）は、もとの1行目の末尾に" : REVERSE COMPLEMENT" とつけ加えたものにし、
かつ、配列部分は 10 文字毎に1文字空白文字を入れ、50 文字毎に改行すること。


＃FASTA フォーマットについては、たとえば以下を参照
＃大文字小文字、曖昧塩基の入力に注意すること
＃　http://fantom2.gsc.riken.go.jp/blast/docs/fasta.html

# 出力はMulti-FASTAとしてパースできればOKです。
# 配列と配列の間に空行を入れるなどしても構いません。

Make a program to output reverse complementary sequence for a multi-FASTA-formatted file.
The 1st line of the output file is represented as original first line with the words ": REVERSE COMPLEMENT".
Add one space for each 10 characters, and insert line feeds for each 50 characters.

＃FASTA format is shown in the following web site.
＃Distinguish capital and small characters. 
＃　http://fantom2.gsc.riken.go.jp/blast/docs/fasta.html
実行例

– 入力ファイルの中身 

>test1
AAaaacGGGGGt
aaaaacGGGGGt
>test2
AAAAAAAGGGwwbbnnnnn

– 出力 

>test1: REVERSE COMPLEMENT
aCCCCCgttt ttaCCCCCgt ttTT
>test2: REVERSE COMPLEMENT
nnnnnvvwwC CCTTTTTTT
▼　３ 配列操作の基礎２（配列の翻訳） 
マルチFASTA形式で複数の塩基配列が書かれたファイルへのパスを引数にとり、それに書かれている塩基配列をアミノ酸配列に翻訳するプログラムを作れ。
"-3"オプションをつけると前向き３フレーム、"-6"で全６フレームの翻訳も行え。
出力はマルチFASTAフォーマットとしてパースできる形式で標準出力へ書き出すこと。

# 出力の配列ヘッダ（＞から始まる行）は実行例にあるように、「元のヘッダ + ": frameN"」のような形式にしてください。
# また、実行例では配列と配列の間に空行はありませんが、空行を入れて見やすくしても構いません。Multi FASTAフォーマットであればOKです。

Make a program which reads sequences in a multi FASTA-formatted file, 
and translates nucleic acid sequences in the file to amino acid sequences.
Path to input file must be accepted as a command line argument.
The program must accept the following two command line option as well: "-3" and "-6".
"-3" and "-6" mean 3-frame translation and 6-frame translation respectively.

# output should be multi FAST-formatted text ,
# and sequence headers (lines starts with ">") should be formatted like "> [original header]: frameN".
実行例
$ cat test.fa
>sequence1
ccg agc aga gct ttc tgg aga gag gaa gag gaa gag gaa gtg gga ggc ggg ccc tag

$ python p3.py test.fa
>sequence1: frame1
PSRAFWREEE EEEVGGGP*

$ python p3.py -3 test.fa
>sequence1: frame1
PSRAFWREEE EEEVGGGP*
>sequence1: frame2
RAELSGERKR KRKWEAGP
>sequence1: frame3
EQSFLERGRG RGSGRRAL

$ python p3.py -6 test.fa
>sequence1: frame1
PSRAFWREEE EEEVGGGP*
>sequence1: frame2
RAELSGERKR KRKWEAGP
>sequence1: frame3
EQSFLERGRG RGSGRRAL
>sequence1: frame -1
LGPASHFLFL FLSPESSAR
>sequence1: frame -2
*GPPPTSSSS SSLQKALL
>sequence1: frame -3
RARLPLPLPL PLSRKLCS

複数のアミノ酸配列が一つのファイルにMulti FASTAフォーマットで収められている。
それらの配列の各々をUniprotのSwissprotに対してblast検索を行い、２番目によく似た配列とのアラインメントのみを取り出すプログラムを作れ。
入力として用いられるMulti FASTAファイルはプログラムの引数として与えられる。
また、blastへは数々のパラメータを与えることができるが、それらは標準設定（何も設定しない）でよい。
出力は次のようなMulti FASTAライクなフォーマットにすること。

（出力例）
--- ここから ----
> sequence1
Query  1    CMVCTV------------SSGPPDRVPGDGSGESCTGYGGVCGASCICGCG---CAVCGC  45 
            ||||||            ||| ||||||           | ||    || |   ||||||     
Sbjct  728  CMVCTVCGDCTGYGASCVSSGRPDRVPG-----------GICG----CGSGESGCAVCGC  772

> sequence2 
Query  1    CIGCRFCSSGF--------ATGQCGCS  19 
            |  || | |||        |   | ||     
Sbjct  81   CESCRHCNSGFLIRNCTVTANAECSCS  107

> sequence3
Query [TAB] [Queryの開始位置] [TAB] [Query配列 (長くても改行を入れない)]
[TAB] [TAB] [QueryとSubjectで塩基が一致している場所には"|"、それ以外は" "]
Subj [TAB] [Subjectの開始位置] [TAB] [Subject配列 (長くても改行を入れない)]
--- ここまで ----

ここで作成されるプログラムは次のように実行できるようにすること。

    $ python p4_noname.py --blast [blastpのパス] --database [DBへのパス] [入力となるMulti-FASTAファイルのパス]

--blast [...]と--database [...]は省略できるようにすること。省略して実行された場合、次のデフォルト値を使うこと。
    --blast: /usr/local/bin/blastp
    --database: ~/blastdb/swissprot

For a given multiple FASTA format sequences, make a command to perform BLAST search to Swissprot database in Uniprot and output only the alignment of the second hit sequence.
参考情報

BLAST は NCBI の ftp サイトからダウンロードして、自分のところにインストールできる．
BLAST+と(Legacy) BLASTの2種類がありますがBLAST+の方を使っておくとよいです。

BLASTの使い方に関しては次のURLなどを参照: http://www.ncbi.nlm.nih.gov/books/NBK1763/

BLAST program can be download from NCBI ftp site and locally installed.

ヒトの全てのRNA配列を含むGenBank形式のファイルを[1]
もしくは[(HGC mirror)]よりダウンロードせよ。
それを用いて、ヒトのmRNAのうち，単純な遺伝暗号による変換では正しいタンパク質配列が得られないものの情報をまとめよ．
以下の処理を行うこと


1.一つのRNA配列についての情報は"LOCUS"で始まる行から "／／" で始まる行までのブロックである．


LOCUS行の最初にIDがかかれており，IDの最初が"NM_"で始まるものがmRNAである．
mRNA配列についてのブロックにのみ処理を行い，ほかは無視すること．



2. "FEATURES"で始まる行以降には，RNAについての様々な情報が記載されている．


そのうち，"CDS"と書かれた行以降には，mRNA配列のうち，蛋白質に翻訳される部分の情報が記載されている．
これらの記述から，各mRNAに対して以下の情報を抽出せよ．

a. CDSの開始および終了の塩基番号（"CDS"と同じ行に記載されている）
b. mRNAが翻訳されたタンパク質のID（"protein_id"=****という形式で記載されている）
c. 翻訳されたタンパク質のアミノ酸配列（translation="*****"という形式で記載されている）



3. "ORIGIN"で始まる行以降にはmRNAの塩基配列が記載されているので，これも抽出せよ
4. 各mRNAのレコードについて，mRNAの塩基配列からCDSで指定された領域を遺伝暗号（問題３参照）に従って翻訳し，タンパク質配列と等しいかどうかを判定せよ.


もし配列が異なる場合にはmRNAのIDとタンパク質のIDをタブ区切りで並べて標準出力へ出力せよ．



5. 4.でみつかったmRNAについて，genbankファイルから該当のブロックを読んで，なぜ翻訳結果が異なるのかを調べよ．


"CDS"の記述のうち，"note"と書かれた部分を参考にするとよい．
 そこで分かった原因をまとめて，4. に加えて出力すること



Download a GenBank formatted file which includes all human RNA sequences from [2] and decompress the file.  
Summarize those mRNAs from which the corresponding protein sequences can not be obtained by translation based on the orthodox genetic code.  

Detailed processes for this problem are as follows, 


1. Information of an RNA sequence is described in a block, which starts with a line starting with "LOCUS" and ending with a line "//".  


Its ID is written just after "LOCUS" identifier, and mRNA IDS starts with "NM_".  
Extract the blocks describing mRNA sequences and ignore the rest.  



2. Lines following a line starting with "FEATURES" include various information about the RNA. 


Among them, lines following a line denoted "CDS" describes the coding region (region which is translated to protein sequence) of mRNA. 




Using these lines, extract following information for each mRNA,

a. Base numbers of the start and the end of CDS (written in the same line as "CDS" sub-header).
b. Protein ID (in a format of "proteinid"=****)
c. Amino acid sequence of the translated protein (written in a format of translation="****")



3. Extract the mRNA sequence from lines following a line starting with "ORIGIN". 
4. For each mRNA block, translate the region in mRNA sequence which was annotated as CDS based on the genetic code (Problem 3),


and compare the translated sequence with the protein sequence obtained in 2-c.  

List the mRNA-ID and protein-ID in a tab-separated style for the entries in which the two amino acid sequences are not identical.  



5. Read the GenBank record of the corresponding to the mRNAs found in 4.  



Analyze why two amino acid sequences differ.  

Hint: Read the part with titled "CDS" and "note".  



Summarize the causes of the inconsistency and report them together with the list obtained in 4.

注意
  • GenBankファイル形式については以下を参照すること: http://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html
  • プログラムの引数としてhuman.rna.gbff.gzを解凍したものへのパスを与えるとファイルを処理して結果を標準出力へ書き出すようにプログラムを書く
    • 自分のプログラムでファイルをダウンロードしたり解凍したりしなくてもOK
  • 上記のファイルリンクは無効なので、[3]からいい感じのを持ってくるとよさげ。(15,4,28 kagaya)
    • コメントありがとうございます。以前は全エントリが一つのファイルに書かれていたのですが、最近では複数ファイルに分割されているようです。
    • human.数字.rna.gbff.gzを全てダウンロードし、それらのファイルを全て一気に処理できるようなプログラムを作成してください。
      • python p4_foo.py /path/to/human.1.rna.gbff.gz /path/to/human.2.rna.gbff.gz ... のように実行できると良いです
▼　６ 立体構造データの解析基礎（PDBMLファイルからアミノ酸配列のFASTAファイルを作る） 
ＰＤＢのファイルを読んで、その中に含まれるタンパク質の配列をFASTA形式（鎖が複数ある場合はマルチFASTA形式）で書き出すスクリプトを作れ。
ＰＤＢの中には２種類の配列（SEQRESレコードとATOMレコード）が書いてあり、通常は一致しているが一致していないこともあるので、
オプション、（-sではSEQRESを読んで、-aではATOMを読んで変換）でどちらの配列からも変換ができるようにせよ。

読み込むべきPDBファイルは、ファイルパスがプログラムの引数として指定される。
また、出力は標準出力へ行え。

プログラムの引数が不正な時、つまり、(i) -s, -a両方が指定された時、(ii) -s, -aどちらも指定されなかった時、または、(iii) PDBファイルのパスが指定されなかった時は、
プログラムが終了ステータス0以外で終了するようにせよ。
プログラムの処理が正常に完了した際は終了ステータス0で終了せよ。

Make a program that reads PDB file and outputs protein sequence as multiple FASTA format.
PDB file includes two types of sequences (SEQRES record, ATOM record). 
Usually these two sequences are identical, but sometimes not. 
Thus, prepare program options (-s: read SEQRES sequence, -a: read ATOM sequence) to deal both sequences.

実行例

– 入力例 (Input example) 

pdb2xhe.ent
(http://www.pdbj.org/pdb_nc/pdb2xhe.ent)

– 出力例 (Output example) 

>2XHEA:UNC18
HMSLKSAVKTVLTNSLRSVADGGDWKVLVVDKPALRMISECARMSEILDLGVTVVEDVSK
QRKVLPQFHGVYFIEPTEENLDYVIRDFADRTPTYEAAHLFFLSPVPDALMAKLASAKAV
KYVKTLKEINTLFIPKEHRVFTLNEPHGLVQYYGSRSSSYNIDHLVRRLSTLCTTMNVAP
IVRYSSTSTPGTERMAMQLQKEIDMSVSQGLINAREGKLKSQFLILDRAVDLKSPLVHEL
TYQAAAYDLLNIENDIYSYSTVDAGGREQQRQVVLGEDDDIWLQMRHLHISEVFRKVKSS
FDEFCVSARRLQGLRDSQQGEGGAGALKQMLKDLPQHREQMQKYSLHLDMSNAINMAFSS
TIDSCTKAEQNIVTEEEQDGNKVRDFIGEVASVVVDRRVSTEDKLRCLMLCVLAKNGTSS
HELNNLLDNANIATPSRSAIYNLEMLGATVVADRRGRKPKTMKRIERDMPYVLSRWTPIV
KDLMEYIATGQLDLESYPAVRDGPSVVQPKRASKSVEEDDDGPATSARKRGNWAKNKGNN
RSLPSTPSGVAVSGNGAAGAAESAKPKLFVFINGTVSYNEIRCAYEVSQSSGYEVYIGAH
NIATPAEFVELVSLLDKADQDVQVLTQGQGDGGLVITTGSAQAGLNLAEV

>2XHEB:SYNTAXIN1
MDRLSRLRQMAAENQPAEASDAAGGAEAQIEETSLSAQPEPFMADFFNRVKRIRDNIEDI
EQAIEQVAQLHTESLVAVSKEDRDRLNEKLQDTMARISALGNKIRADLKQIEKENKRAQQ
EGTFEDGTVSTDLRIRQSQHSSLSRKFVKVMTRYNDVQAENKRRYGENVARQCRVVEPSL
SDDAIQKVIEHGTEGIFSGMRLEGAEAKLNEIRDRHKDIQQLERSLLELHEMFTDMSTLV
ASQGEMIDRIEFSVEQSHNYVKKATEQVVQARHYQESAR
▼　７ 立体構造データの解析応用（PDBMLファイルとPDBの座標の扱い） 
PDBMLファイルを読んで重心の座標を表示するプログラムを作れ。
なお、座標ファイルはサイズが大きくなるため通常gzip形式で圧縮してある。
作成するスクリプトでは、gzip圧縮されたファイルと無圧縮のファイルをそのまま扱えるようにせよ。
出力は次のような形式で標準出力へ行え。

---- ここから ----
x [タブ] [重心のx座標] [改行]
y [タブ] [重心のy座標] [改行]
z [タブ] [重心のz座標] [改行]
---- ここまで ----

# gzip圧縮されたファイルのファイル名は拡張子".gz"で終わる。それ以外の拡張子であれば無圧縮ファイルと判断してよい。

注意

問題５で用いたＰＤＢファイルではなく、ＸＭＬ形式のファイルを入力とすること。

PDBMLでは原子ひとつが <PDBx:atom_site> タグひとつに相当する。 

        <PDBx:atom_site id="1"> ... </PDBx:atom_site>

原子座標については下記タグを見ればよい。 

         <PDBx:Cartn_x> </PDBx:Cartn_x>
         <PDBx:Cartn_y> </PDBx:Cartn_y>
         <PDBx:Cartn_z> </PDBx:Cartn_z>

PDBMLファイルの読み込みは通常のXMLパーサー（xml.saxなど）を使用すると良い。

PDBMLファイルの例を下に示す。 

<?xml version="1.0" encoding="UTF-8" ?>

<PDBx:datablock datablockName="1ZCB"
   xmlns:PDBx="http://pdbml.pdb.org/schema/pdbx-v32.xsd"
   xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
   xsi:schemaLocation="http://pdbml.pdb.org/schema/pdbx-v32.xsd pdbx-v32.xsd">

   <PDBx:atom_siteCategory>
      <PDBx:atom_site id="1">
         <PDBx:B_iso_or_equiv>61.18</PDBx:B_iso_or_equiv>
         <PDBx:B_iso_or_equiv_esd xsi:nil="true" />
         <PDBx:Cartn_x>-55.369</PDBx:Cartn_x>
         <PDBx:Cartn_x_esd xsi:nil="true" />
         <PDBx:Cartn_y>-2.728</PDBx:Cartn_y>
         <PDBx:Cartn_y_esd xsi:nil="true" />
         <PDBx:Cartn_z>36.394</PDBx:Cartn_z>
         <PDBx:Cartn_z_esd xsi:nil="true" />
         <PDBx:auth_asym_id>A</PDBx:auth_asym_id>
         <PDBx:auth_atom_id>N</PDBx:auth_atom_id>
         <PDBx:auth_comp_id>ALA</PDBx:auth_comp_id>
         <PDBx:auth_seq_id>46</PDBx:auth_seq_id>
         <PDBx:group_PDB>ATOM</PDBx:group_PDB>
         <PDBx:label_alt_id></PDBx:label_alt_id>
         <PDBx:label_asym_id>A</PDBx:label_asym_id>
         <PDBx:label_atom_id>N</PDBx:label_atom_id>
         <PDBx:label_comp_id>ALA</PDBx:label_comp_id>
         <PDBx:label_entity_id>1</PDBx:label_entity_id>
         <PDBx:label_seq_id>31</PDBx:label_seq_id>
         <PDBx:occupancy>1.00</PDBx:occupancy>
         <PDBx:occupancy_esd xsi:nil="true" />
         <PDBx:pdbx_PDB_ins_code xsi:nil="true" />
         <PDBx:pdbx_PDB_model_num>1</PDBx:pdbx_PDB_model_num>
         <PDBx:pdbx_formal_charge xsi:nil="true" />
         <PDBx:type_symbol>N</PDBx:type_symbol>
      </PDBx:atom_site>
      <PDBx:atom_site id="2">
         <PDBx:B_iso_or_equiv>61.01</PDBx:B_iso_or_equiv>
         <PDBx:B_iso_or_equiv_esd xsi:nil="true" />
         <PDBx:Cartn_x>-56.584</PDBx:Cartn_x>
         <PDBx:Cartn_x_esd xsi:nil="true" />
         <PDBx:Cartn_y>-2.646</PDBx:Cartn_y>
         <PDBx:Cartn_y_esd xsi:nil="true" />
         <PDBx:Cartn_z>35.524</PDBx:Cartn_z>
         <PDBx:Cartn_z_esd xsi:nil="true" />
         <PDBx:auth_asym_id>A</PDBx:auth_asym_id>
         <PDBx:auth_atom_id>CA</PDBx:auth_atom_id>
         <PDBx:auth_comp_id>ALA</PDBx:auth_comp_id>
         <PDBx:auth_seq_id>46</PDBx:auth_seq_id>
         <PDBx:group_PDB>ATOM</PDBx:group_PDB>
         <PDBx:label_alt_id></PDBx:label_alt_id>
         <PDBx:label_asym_id>A</PDBx:label_asym_id>
         <PDBx:label_atom_id>CA</PDBx:label_atom_id>
         <PDBx:label_comp_id>ALA</PDBx:label_comp_id>
         <PDBx:label_entity_id>1</PDBx:label_entity_id>
         <PDBx:label_seq_id>31</PDBx:label_seq_id>
         <PDBx:occupancy>1.00</PDBx:occupancy>
         <PDBx:occupancy_esd xsi:nil="true" />
         <PDBx:pdbx_PDB_ins_code xsi:nil="true" />
         <PDBx:pdbx_PDB_model_num>1</PDBx:pdbx_PDB_model_num>
         <PDBx:pdbx_formal_charge xsi:nil="true" />
         <PDBx:type_symbol>C</PDBx:type_symbol>
      </PDBx:atom_site>
      <PDBx:atom_site id="3">
         <PDBx:B_iso_or_equiv>60.52</PDBx:B_iso_or_equiv>
         <PDBx:B_iso_or_equiv_esd xsi:nil="true" />
         <PDBx:Cartn_x>-56.371</PDBx:Cartn_x>
         <PDBx:Cartn_x_esd xsi:nil="true" />

応用的課題

引数として与えられた文字列パターンを含むエントリーをUniprotのSwissprotデータベースから検索し、マッチのあったエントリー名とパターンの１文字目の位置を出力せよ。
これを使って、”LSLAVAG”という部分配列がいくつあるか調べよ。また、任意の文字とマッチする文字Ｘを指定できるように拡張せよ。
Uniprotはgw.hgc.jp:/usr/local/db/flat/uniprot/以下にある。

DBM の機能を使って、Uniprot のSwissprotのエントリー名を与えると、即座にその内容を出力するコマンド
sgetを作れ（別のコマンドを作って、あらかじめ各エントリーの開始位置をデータベース化しておく）。

GenBank Primate
divisionについて、それぞれデータベースのエントリー数、核酸配列長の平均値、標準偏差、合計値、最大値、最小値を求めよ。最大最小値については、エントリー名とその定義行も出力せよ。また、'*'
文字を頻度に応じて並べた長さの分布図を作れ。

# GenBank については、
# /usr/local/db/flat/genbankの下にある
# Primate は gbpri*.seq の形のもの
# 分布図はたとえば以下のような形で（別に対数スケールをとる必要はないが）
# Genbankの処理にはBioperlを使うと良い。

frequency 1   10   100 1000 10000
length(aa) ++----+----+----+----+----+
 0- 500 |*************************  80018
500-1000 |**********************  17541
1000-1500 |******************  2782
1500-2000 |***************  653
2000-2500 |*************  327


#以下の手順で遺伝子情報(EntrezGeneD，RefseqID, Symbol, コーディング領域の遺伝子座標)がまとまったビューを作成するコマンドを作れ．
 - refFlat.txt, gene2refseq, geneinfoの各データからSQLを使用し，テーブルを作成せよ．
 - SQLにより，3種のテーブルを結合し，ビューを作成せよ．

2.1で作成したビューに対して，条件(染色体番号，開始点，終了点)を満たす遺伝子の情報を返すコマンドを作れ．
- 染色体番号は一致するものを取得する．
- コーディング領域に関しては，検索条件で指定される遺伝子座標の範囲内のものを全て取得する．


# pythonのsqlite3モジュールを使用して，データベースに対する操作をSQLを行う．
# 遺伝子座標は，染色体番号と開始点・終了点のことを指す．
# 各種データは膨大なので，Homo_sapiensのものだけで良いと思います

– 検索結果例 染色体番号:4，開始点:600000，終了点:700000 

# refseq_id, entrez_gene_id, symbol, chrom, cds_start, cds_end
NM_000283	5158	PDE6B	chr4	619415	663896
NM_001145291	5158	PDE6B	chr4	619415	663896
NM_001145292	5158	PDE6B	chr4	647766	663896
NM_001294341	84179	MFSD7	chr4	675746	682916
NM_001294342	84179	MFSD7	chr4	675746	682916
NM_002477	4636	MYL5	chr4	671815	675783
NM_007100	521	ATP5I	chr4	666288	668036
NM_032219	84179	MFSD7	chr4	675746	682916
NR_033743	521	ATP5I	chr4	668127	668127

応用的課題 (理論編)

(バイオインフォマティクス -ゲノム配列から機能解析へ- 第2版 p.111 第3部 配列アラインメントスコアを計算する)

BLOSUMの方法による対数オッズとオッズスコアの計算 

類縁のあるよく似た配列セットのアラインメントの1つの列において、アミノ酸Dはアミノ酸Eに頻度0.10で置換し、
この列においてDとEの頻度に基づくこの置換が期待される回数は0.05である。

#アラインメントにおけるDからEへの置換が見つかるオッズスコアはいくつか。

2. ビット単位でのDからEへの置換の対数オッズスコアはいくつか。

3. この置換に相当するBLOSUMアミノ酸スコア行列中のエントリは何か。
　結果をBLOSUM62行列の実際のエントリと比較せよ。

4. 同じ列においてDが変異しない割合は0.80で、Dが変異しない割合の期待値は0.10であった。
　相当する対数オッズスコアとDが変異しないBLOSUM62のエントリを計算せよ。

短いアラインメントに対する対数オッズとオッズスコア 

#上記の値を用いて以下のアラインメントの対数オッズスコアはビット単位でいくつになるか。

      DEDEDEDE
      DDDDDDDD

2. 上のアラインメントのオッズスコアはいくつか。

(バイオインフォマティクス -ゲノム配列から機能解析へ- 第2版 p.112 第4部 小さいもしくは大きいギャップペナルティでアラインメントスコアを計算する) 

次のような操作を行い、ギャップペナルティとアラインメントスコアの関係について述べよ。

#UniProtから以下の配列を得る。
　* E. coliのrecA.pro (P0A7G6) (テキストではP03017だが、obsoleteになっているので変更)
　* Saccharomyces cerevisiaeのrad51.pro (P25454)
　　このタンパク質は同じ機能を持っている。それは相同な一本鎖DNAが対合するのを促進する機能である。
　　それらは確実にほぼ同じ立体構造を持っているが、アラインメントするのが難しいくらい十分分岐している。

2. LALIGNを次のようなギャップペナルティの設定で実行し、得られたアラインメントの特徴を記録する。
**ギャップペナルティ(open, ext) = -12, -2
**ギャップペナルティ(open, ext) = -5, -1

講義

Copyright © 木下・西研究室｜東北大学大学院・情報科学研究科生命情報システム科学分野 All Rights Reserved.