M1問題/M1 exercise

 

以下の問題を解き、報告会で解答を報告してください。使用言語は Python がおすすめです。
Please doing the following exercises and report the answers at the meeting. Recommended language is Python.

1. アミノ酸配列からRNA 配列への変換/Transform an Amino Acid Sequence into RNA sequences

問題/Problem: Homeobox の antennapedia 型タンパク質はそのアミノ酸配列中に LFPWM という配列パターンを持つことが知られている。アミノ酸配列 LFPWM をコードする RNA 配列のパターンを全て列挙せよ。遺伝コードは Universal genetic code に従うものとする。なお、 LFPWM に限らずあらゆるアミノ酸配列に対応できるようにプログラムを作成せよ。
Homeobox’s antennapedia type protein has a amino acid sequence pattern “LFPWM”. List all possible RNA sequence patterns coding this amino acid sequence pattern. Genetic code shall conform to Universal genetic code. You have to make a program suitable for not only “LFPWM” but also all possible amino acid sequences.

2. 制限酵素サイトの検索/Search Restriction Enzyme Sites

準備/Preparation: λファージの完全ゲノム配列をテキストファイルに保存してください。次の RefSeq ページ(NC_001416)で FASTA をクリックすると FASTA 形式のゲノム配列を表示できます。
Please save the whole genome sequence of λ phage in a text file. Click “FASTA” on this RefSeq page (NC_001416) and get the genome in FASTA format.

問題/Problem: λファージの完全ゲノムを制限酵素 EcoRI で完全に消化したとき、得られる全ての断片の長さを答えよ。HindIII, BamHI, NotI で消化した場合についても同様に答えよ。制限酵素サイトは REBASE で検索せよ。
When digesting the whole genome sequence of λ phage completely with EcoRI, how long can each fragment be obtained? How about HindIII, BamHI and NotI? For restriction enzyme sites, please refer to REBASE.

3. データベースの現状/Current state of database (An access to the supercomputer is required)

問題/Problem 3-1:

1) KEGG GENOME の塩基配列データ(スパコン:/bio/db/fasta/genome/)と KEGG GENES の遺伝子配列データ(スパコン:/bio/db/fasta/genes/)を読み込み、各ゲノム(T番号毎)の総塩基数とタンパク質コーディング遺伝子の数をカウントせよ。また、横軸をゲノムサイズ、縦軸をタンパク質コーディング遺伝子の数とするプロットを作成せよ。
Read the base sequence data of KEGG GENOME (supercomputer:/bio/db/fasta/genome/) and the gene sequence data of KEGG GENES (supercomputer:/bio/db/fasta/genes/) and count the total number of bases and protein-coding genes of each genome (T number). Also, create plots with the genome size on the horizontal axis and the number of protein-coding genes on the vertical axis.

2) 次に、 スパコン:/bio/db/ideas/genome/ から生物種名、ゲノムサイズ、遺伝子数を取得して表を作成し、上記のカウントと合うかを確かめよ。
Next, obtaine species name, genome size, number of genes from supercomputer:/bio/db/ideas/genome/ and create the table of them. Also, check if the genome size and number of genes of each genome matches the above count.

問題/Problem 3-2: スパコン:/bio/db/ideas/ にはゲノムネットでサービスしているデータベースがある。このうち、以下のデータベースについて、そのエントリー数をカウントして、表にせよ。また、それぞれのデータベースの内容について簡単にまとめよ。
supercomputer:/bio/db/ideas/ has a database serving on the GenomeNet. Among them, count the number of entries for the following databases and tabulate them. Also, report what kind of contents databases have. (eg. “RefSeq is manually curated database for non-redundant set of sequences provided by NCBI. That includes genomic DNA, transcripts, and proteins.” Please report like this!!)

RefSeq, PDB, SWISS-PROT, TrEMBL, PFAM, NCBI-CDD, PATHWAY, MODULE, KO, COMPOUND, ENZYME, REACTION, DRUG, DISEASE

4. EC 番号のソート/Sort EC number (An access to the supercomputer is required)

問題/Problem: KEGG KO(スパコン:/bio/db/ideas/ko/)をパースし、DEFINITION 行に EC 番号があるものを抽出し、抽出した EC 番号を、各カラムの数字の大きさに従ってソートして、EC 番号と対応する KO ID の一覧表として作成せよ。
Parsing the latest version of KEGG KO (supercomputer:/bio/db/ideas/ko/), extracting those with an EC number in the DEFINITION line, sorting the extracted EC numbers according to the size of the numbers in each column, and creating them as a list of EC numbers and corresponding KO IDs.

5. DBGET/BLASTデータベースの作成/Create DBGET/BLAST database (An access to the supercomputer is required)

問題/Problem 5-1: DBGET の bfind コマンドを用いて KEGG GENES(スパコン:/bio/db/fasta/genes/)から自分の興味のある遺伝子(例:ribosome)をキーワード検索で抜き出せ。そして、seqnew コマンドを用いて DBGET データベースを作成せよ。作成したデータベースから bget コマンドでエントリを取得できることを確認せよ。
Extract the gene of your interest (eg. ribosome) from KEGG GENES (supercomputer:/bio/db/fasta/genes/) by keyword search using bfind command in DBGET. Then, create DBGET database using seqnew command.
From the database you created, make sure you can get entries with bget command.

問題/Problem 5-2:

1) DBGET の bfind コマンドを用いて KEGG GENES から自分の興味のある遺伝子をキーワード検索で抜き出し、FASTA 形式で保存せよ。そして、BLAST+ の makeblastdb を用いて BLAST データベースを作成せよ。作成したデータベースから blastdbcmd コマンドでエントリを取得できることを確認せよ。
Extract the gene of your interest from KEGG GENES by keyword search using seqnew command and save it in FASTA format. Then, create BLAST database using makeblastdb command in BLAST+. From the database you created, make sure you can get entries with blastdbcmd command.

2) 次に、好きな配列をクエリとして、作成したデータベースに対して BLAST 検索し、上位10個の配列のFASTAファイルを作成せよ。
Next, Using the sequence of your interest as a query, search BLAST against the created database and make the FASTA file of top 10 sequences.

[DBGET help] -h オプションで簡略な、man コマンドで詳細なヘルプを表示できるので活用してください。
Please see help as it appears in -h option (brief help) or man command (detail help).

[BLAST+ help] -h オプションで簡略な、-help コマンドで詳細なヘルプを表示できるので活用してください。下の公式マニュアルも役に立ちます。
Please see help as it appears in -h option (brief help) or -help option (detail help). The official manual below is also useful.
BLAST command line applications user manual (pdf)

6. BLASTの結果ファイルのパース/Parsing the BLAST result file

準備/Preparation: 大腸菌 O157:H7 堺株のリボソーム生合成タンパク質 YsxC のアミノ酸配列をテキストファイルに保存してください。次の KEGG GENES ページ(ecs:ECs4787)で AA seq をクリックすると FASTA 形式のアミノ酸配列を表示できます。大腸菌 K-12 MG1655 の転移 RNA リシジン合成酵素のアミノ酸配列(eco:b0188)も同様にテキストファイルに保存してください。
Please save in a text file the amino acid sequence of ribosome biogenesis protein YsxC of E. coli O157:H7 str. sakai. Click “AA seq” on this KEGG GENES page (ecs:ECs4787) and get the amino acid sequence in FASTA format. Please save also the amino acid sequence of tRNA-lysidine synthetase of E. coli K-12 MG1655 (eco:b0188).

問題/Problem: リボソーム生合成タンパク質 YsxC のアミノ酸配列をクエリーとして、 KEGG GENES データベースに対して BLAST 検索を行い、結果を TSV 形式で保存せよ( E-value の閾値は 1e-10 とせよ)。保存したTSV 形式をHTML 形式の一覧表として表示せよ。転移 RNA リシジン合成酵素についても同様にせよ。(BLAST コマンドの HTML 出力オプションは使わないでください)
Search the KEGG GENES database with BLAST for the amino acid sequence of ribosome biogenesis protein YsxC as a query and save the result as TSV format (E-value threshold should be 1e-10). Then, display the result file in HTML format. Do the same with tRNA-lysidine synthetase. (Please do not use a HTML output option of BLAST commands)

7. PATHWAY データからのネットワーク特徴量計算/Calculation of network feature quantity from PATHWAY data

準備/Preparation: 大腸菌 の KEGG PATHWAY データを下記リンクからダウンロードしてください(研究室のLANからのみアクセス可能)。
Please download the KEGG PATHWAY dataset of E. coli from the link below (available only from lab’s LAN).
KEGG PATHWAY xml files of E. coli (tar ball)

問題/Problem:

1) 大腸菌 KEGG PATHWAY データ(XML ファイル)から化合物の二項関係(反応の基質と生成物の関係)を抽出せよ。
Extract the binomial relationship of the compound (the relationship between the substrate of the reaction and the product) from the KEGG PATHWAY dataset of E. coli (XML files).

2) 次に、各化合物の次数(いくつの化合物と関係があるか)を計算し、そのヒストグラムを作成せよ。
Next, calculate the degree of each compound (how many compounds are related) and create a histogram of that.

3) そして、抽出した二項関係が構成するネットワークにいくつの連結成分があるかを計算せよ。また、各連結成分のサイズ(何個の化合物が含まれるか)も計算せよ。その結果から、サイズごとの連結成分の個数をヒストグラムとして表せ。
Then, calculate how many connected components are in the network constructed with binary relation. Also calculate the size of each connected component (how many compounds in it). From the result, create a histogram of the size of the connected component.

8. 16S rRNA 配列解析/16S rRNA sequence analysis

準備/Preparation: 下記リンクにアクセスし(パスワードは virus です)、ディレクトリ Data/14Aug/16Sread_merged に移動してください。そして、ディレクトリ内の9個の FASTA ファイルをダウンロードしてください。
Please access the link below (password is “OBV”) and move to the directory Data/14Aug/16Sread_merged. Then, please download the nine FASTA files there.
Osaka Bay 16S merged reads

問題/Problem:

1) 16S rRNA 配列の FASTA 形式データ(9サンプル)を読み込んで、マージして一つのファイルにせよ。そして、CD-HIT を使用して 97% identity でクラスタリング(OTU 化)せよ。ここで使用する9つのサンプルは大阪湾から3時間おきに取水したものである。B1 は12:00、B2 は15:00、B3 は18:00・・・そして B9 は次の日の12:00の取水である。
Read FASTA format data (9 samples) of 16S rRNA sequence and merge them into one file. Then, cluster it (turn to OTU) with 97% identity using CD-HIT. The nine samples used here are taken from Osaka bay every three hours. B1 was taken at 12:00, B2 at 15:00, B3 at 18:00… and B9 at 12:00 of the next day.

2) CD-HIT を使用して、SILVA(/bio/db/fasta/ribosomaldb/)をリファレンスとして各 OTU の代表配列に生物分類をアノテーションせよ。ここでも、97% identity 以上を閾値とする。
Annotate biological classification for representative sequences of each OTU. Srarch against SILVA (/bio/db/fasta/ribosomaldb/) reference database using CD-HIT with 97% identity.

3) アノテーションした結果から OTU の各サンプルにおける規格化された相対(正規化)頻度の表を作成せよ。そして、Phylum および Class レベルでのサンプル毎の分類学的組成を示す積み上げ棒グラフを作成せよ。
From the annotated results, make the the table of relative (normalized) frequency of OTUs in each sample. Then, create a stacked bar graph that shows the taxonomic composition of each sample in Phylum and Class level.

4) 相対頻度の表から、Bray-Curtis Dissimilarity を距離関数として用いて、サンプル間の距離行列を計算せよ(OTUレベルで)。そして、その距離行列に多次元尺度構成法(MDS)を適用して、2次元および3次元空間にプロットせよ。Bray-Curtis Dissimilarity と MDS の計算には、それぞれ R の vegan パッケージと stats パッケージを使うと便利である。
From the table of relative frequency, calculate a distance matrix between samples using Bray-Curtis Dissimilarity as distance function (in OTU level). Then apply multidimensional scaling method (MDS) on the distance matrix and plot samples on two- and three-dimensional space. For calculation of Bray-Curtis Dissimilarity and MDS, it is convenient to use R vegan package and stats package respectively.

[CD-HIT help] CD-HIT は様々な用途に対応しており、実行プログラム名 (cd-hit-xxx) は用途によって異なります。下記ユーザーガイドをよく読んで、適切なプログラムを選択しましょう。
CD-HIT supports various applications, and the execution program name (cd-hit-xxx) varies
depending on the application. Please read the User’s Guide below carefully and select the appropriate program.
CD-HIT user’s guide

[Reference] Yoshida T, et al. Locality and diel cycling of viral production revealed by a 24 h time course cross-omics analysis in a coastal region of Japan. 2018. Isme J 12(5):1287-95 doi:10.1038/s41396-018-0052-x

9. Mollivirusのトランスクリプトーム解析/Transcriptomic analysis of Molliivirus

準備/Preparation: NCBIの Sequence Read Archive (SRA) で SRX1078581 と検索すると、この問題で使用するMollivirusのRNA-seqデータが3つ表示されます。これらのデータを SRA Toolkit の prefetch コマンドを使用してダウンロードし、fastq-dump コマンドで FASTQ 形式に変換してください。またMollivirusの参照ゲノム(NC_027867)並びに宿主 Acanthamoeba castellanii のゲノムアセンブリ(GCF_000313135.1)を FASTA 形式でダウンロードし、これら2つのゲノムを1つのファイル(multi-FASTA 形式)に統合してください。

You can find RNA-seq data of Mollivirus with ID:SRX1078581 in NCBI Sequence Read Archive (SRA). Then please get those data using prefetch command of SRA Toolkit and convert them into FASTQ format with fastq-dump command. Next, please download the referance genome of Mollivirus (NC_027867) and a genome assembly of its host Acanthamoeba castellanii (GCF_000313135.1) in FASTA format and merge these two genomes in one file (i.e. multi-FASTA format).

問題/Problem:

1) それぞれのシークエンスデータをMollivirusとアメーバの統合リファレンスゲノムにマッピングせよ(ツール例:HISAT2、STAR)。マップされたリード数、マップされたリードの割合、Mollivirusゲノムにマップされたリード数、Mollivirusにマップされたリード数とアメーバゲノムにマップされたリード数の比率を算出してグラフで表示せよ。
Do genome mapping against integrated reference genome of Mollivirus and A.castellanii. calculate the number of mapped reads, the ratio of mapped reads, the number of reads mapped to Mollivirus genome, the number of reads mapped to A.castellanii and the ratio of reads mapped to Mollivirus or A.castellanii.

2) Integrative Genomics Viewer (IGV) でマッピング結果を可視化できる。発現していると思われるMollivirusの遺伝子を一つ選び、遺伝子の発現量を時系列順に並べて表示し、そのスクリーンショットを示せ。スパコンでIGVを起動するにはsshコマンドで-Yオプションが必要になる。
You can visualize the result of genome mapping. please choose a annotated gene and show the results in the order of time scale. You can use IGV on the supercomputer, but you need -Y option with ssh command.

3) アノテーションファイルを使用して、2)で選んだ遺伝子の転写量をFPKMで正規化してグラフ表示せよ。
Calculate FPKM of the gene you chose.

[Reference] Legendre M, et al. In-depth study of Mollivirus sibericum, a new 30,000-y-old giant virus infecting Acanthamoeba. Proc Natl Acad Sci U S A. 2015. 112(38):E5327-35 doi:10.1073/pnas.1510795112

10. ワインの質の機械学習による予測/Predict quality of wine with machine learning

準備/Preparation: Wine Quality Datasetsよりwinequality.zipをダウンロードして解凍してください。この問題ではデータセットwinequality-white.csvを使います。このデータセットはポルトガルのワイン「ヴィーニョ・ウェルデ」(Vinho Verde)の白ワインボトルの評価の結果です。第1列から第11列までは物理化学量の測定値、第12列は専門家によるワインの質の評価です。この問題を解くには、Python ユーザなら機械学習パッケージ scikit-learn を使うと便利です。
Please download “winequality.zip” from Wine Quality Datasets and unzip it. A dataset “winequality-white.csv” is used in this problem. This dataset is the result of evaluation of bottles of Portuguese white wine “Vinho Verde”. The 1st column to the 11th columns are meseurements of the physicochemical parameters and the 12th column is the quality of wine evalueted by experts. The machine learning package “scikit-learn” is convinient for solving this problem if you are Python user.

問題/Problem:

1) 白ワインデータセットの散布図行列を作れ。散布図行列により特徴量同士の相関関係をまとめて可視化できる。
Create a scatterplot matrix of the white wine dataset. The scatterplot allows us to visualize the pair-wise correlation between the different features in one place.

2) ワインの質を物理化学量から予測したい。まずワインの質の列を、7から10までを良いワイン(陽性クラス)、1から6までを悪いワイン(陰性クラス)の2クラスに変換せよ。次に、全サンプルの70%を学習セット、30%をテストセットに分割せよ。そして、物理化学量の列を入力特徴量として用いてワインの質(良いか悪いか)を予測するサポートベクトルマシン(SVM)を学習せよ。テストセットでの正答率と混同行列を報告せよ。入力特徴量は標準化すること(平均を0、標準偏差を1に合わせる)。ハイパーパラメータはグリッドサーチによりチューニングすること。
We want to predict the quality of wine based on the physicochemical parameters. First, convert the column of quality to two classes like: 7 to 10 is good wine (a positive class) and 1 to 6 is bad wine (a negative class). Next, split the dataset into a training set (70% of all samples) and a test set (30% of all). Then, use columns of the physicochemical parameters as input features to train a Support Vector Machine (SVM) to predict the quality of wines (good or bad). Report accuracy and a confusion matrix on the test set. Input features should be standardized (centered at mean zero with standard deviation one). Hyperparameters should be tuned via grid search.

3) テストセットでの受信者動作特性(ROC)曲線を描け。ROC曲線の計算には決定関数の返り値を使う。
Plot a ROC (Reciever Operating Characteristic) curve on the test set. Return values of decision function will be used for creating a ROC curve.

4) 各入力特徴量の相対重要度を計算し、棒グラフにせよ。相対重要度は次のように定義される。
Calculate relative importance of each input feature and plot bar graph of them. Relative importance is defined as below.
相対重要度の定義/Definition of relative importance: Definition of relative I
importance (pdf)

[Reference] Cortez P, et al. Modeling wine preferences by mining from physicochemical properties. 2009. Decis Support Syst 47:547-53 doi:10.1016/j.dss.2009.05.016