未回答質問

アミノ酸配列でサイトごとに計算される保存度は何を意味しているのか？

study_bioinfo — Mon, 30 Sep 2019 04:10:08 +0900

参考書に載っていた内容

・どうして保存度を算出するのか？

→アラインメントから構造や機能の情報を抽出するため

・サイトごとに計測する理由は？

→サイトレベルで情報を抽出する方法により、機能あるいは構造と関連があると推測されるサイトについては、アミノ酸置換などの突然変異導入実験により、その予測の確認を比較的容易に行うことができる。モチーフの同定。

・構造とは？

→立体構造のこと

・機能とは？

→タンパク質が折りたたまれ、他のタンパク質等と相互作用を起こすこと。

そのため、立体構造が似ていることは同様の機能を有していることを予測することができる。

マルチプルアラインメントされたアミノ酸配列で、サイトごとに保存度を「以下の式」で計算する意義について知りたいです。

保存度は以下の式で計算されるものとします。

・Valder氏の保存度算出式(https://gyazo.com/8285e51d7bd7db8d8057bb376e9e454f)

・スコアマトリックスにはBlosum62を使用(https://gyazo.com/7c8c0bb28db10ba26d07bdf6e40aed57)

・スコアマトリックスの変換式に、写真の変換式を使用(https://gyazo.com/60242cf27db66929a5490dc0ed9b1558)

・重み付けにはHenikoff氏の重み付けを使用(https://gyazo.com/7907a5beddc2a1ee6348e310f9a5d78b)

保存度を計算する理由について、私の考察

・完全保存されているサイト(=モチーフ)を発見すること。

理由：Valder氏の考案した式によれば、完全保存されているサイトについては、どの配列で比較しても、一意の値に収束するから。反対に、完全保存されていなければ、一意の値に収束することはなく、値にばらつきが見られるから。

考察に対する反論

完全保存されているサイトを発見するためだけであれば、別に保存度が必要ではない。視覚的に完全に保存されているサイトを見つければいい。

よって、私の考察とは異なる目的があるのだろう。という結論にたどり着きました。ご回答いただけると幸いです。

pyBedtoolsによるベン図の作成

Yoda — Thu, 01 Nov 2018 06:53:15 +0900

またお願いします、DRY解析教本にてChip seqの解析を勉強している者です。

venn_mpl.py~コマンドを使用してベン図を作成しようとしておりますが、”matplotlib is required to make a Venn diagram with venn_mpl.py”と表示され解析できません。 Homebrewでmatplotlibをインストールしても同じ結果です。

お知恵をお貸し下さい、よろしくお願いいたします。

BAMファイルのソートについて

Yoda — Wed, 31 Oct 2018 07:57:48 +0900

度々失礼します、次世代シークエンサーＤＲＹ解析教本でCHIPseq解析の勉強をしている初心者です。

bowtieでマッピング後、samファイルをbamファイルに変換しました。次に教科書通り「samtools sort 検体.bam 検体_sorted」というコマンドでソートをしようとしてもできません。そこでネットを参考に「samtools sort -@ 4 検体.bam 検体_sorted」でやってみましたが、検体_sortedファイルは作成されるのですが肝心の検体_sorted.bamファイルは作成されず、次のindex作成に進めず困っています。

何度もすいません、ご教授よろしくお願いいたします。

ngsplotdbのインストールについて

Yoda — Sat, 20 Oct 2018 02:58:09 +0900

次世代シークエンサーＤＲＹ解析教本でCHIPseq解析の勉強をしている初心者です。パソコンはMACBOOK PRO（macOS High Sierra）を使用しています。

教科書通りGoogle Drivehttps://drive.google.com/drive/folders/0B1PVLadG_dCKNEsybkh5TE9XZ1Eよりngsplotdb_mm9_67_3.00.tar.gzとngsplot_mm9_67_3.00_enhancer.tar.gzをダウンロードしngsplotdb.py install〜というコマンドでインストールをしようとしています。

まずパソコンにダウンロードした時点でどうしてもtar.gzファイルがtarファイルに変換（？）されてしまい、gzが抜けてしまいます。その後ngsplot.py installコマンドでインストールしようとしても、Downloaded file may be corruptedと表示されインストールすることができません。

どなたか解決策がわかるかたお願いいたします。

DESeq2でtranscripts IDとGene IDを紐付けする方法

mikan03 — Sat, 28 Jul 2018 17:05:55 +0900

いつもお世話になっています。

おかげさまでhisat2 -> stiringTie -> DESeq2まで進みそれらしいデータも取れました。

DESeq2ではTranscript_IDでデータが出てきますが、これをGene_IDに紐付けすることができません。 Mergeなど試してみましたが、data0(O obs)が返されます。解決策をご教示いただければ幸いです。

ーーーーー

たとえば部分的なファイルを切り出してmergeしてみるとうまく行きます。

trancript_ID33 <- c("NR_131893","NR_131893")

gene_name33 <- c("Porcn", "Porcn")

gene_id33 <- c("MSTRG.71262","MSTRG.71262")

test33 <- data.frame(trancript_ID33,gene_name33,gene_id33) 　 baseMean <- c(53.33436)

log2FoldChange <- c(-0.431138)

pValue <- c(0.2859373)

gene_id_NR <- c("NR_131893")

test66 <- data.frame(baseMean,log2FoldChange,pValue,gene_id_NR)

テスト材料作ってからmerge

merge99 <- merge(test33,test66)

これはうまくできました。

全体でMergeするとうまくいかないところなのですが、なんとかする方法はありますでしょうか？バイオではなくてRの問題かもしれませんが、ご教示いただけますと幸いです。よろしくお願いします。（そもそも根本的に間違っているかもしれません。）

追記：Rでの作業を書きます。

library("DESeq2")

公式ページに従って入力。transcritp_count_matrix.csvはstringtieで生成

countData <- as.matrix(read.csv("transcript_count_matrix.csv", row.names="transcript_id"))

合計6未満(N=3で２郡間の比較）の微弱発現を除去

countData1 <- countData[apply(countData,1,sum)>6,]

PHENO_DATA.txtでサンプル割り当て。#PHENO_DATA.txtの中身

Genotype

1 WT

2 WT

3 WT

4 KO

5 KO

6 KO

colData <- read.csv("PHENO_DATA.txt", sep="t", row.names=1)

確認1

all(rownames(colData) %in% colnames(countData1))

[1] TRUE がでるのでOK

all(rownames(colData) == colnames(countData1))

[1] TRUE がでるのでOK

DEseq2に流し込み。

dds <- DESeqDataSetFromMatrix(countData = countData1, colData = colData, design = ~ Genotype) dds <- DESeq(dds) res <- results(dds) res_gene_id$gene_id <- row.names(res)

Gene IDをつけるためrtracklayarうごかす。

Stringtieでmergeしたファイルをいれる。

library(rtracklayer) gtf <- readGFF("stringtie.merge.gtf")

stringtieの場合はtypeがexonかtranscript のようなので、exonにする。

gtf_gene <- subset(gtf, gtf$type == "exon")

紐付けされたリストを取り出しておく。

gtf_gene_1 <- gtf[,c("gene_id","gene_name","transcript_id")] colnames(gtf)

**#二つの票を合体することで、transcript_IDとgeneIDを紐付けする。ここがうまくいかない

DEseq2_result_ref <- merge(DEseq2_result,gtf_gene_1,all=T, sort=F)**

CSVへの書き出しで終わるはず。

write.csv(DEseq2_result_ref,"DEseq2_result_ref.csv", quote=FALSE, row.names=FALSE)

参考ページ(どうもありがとうございます。）

https://ncrna.jp/blog/item/388-deseq2-ggplot2

http://ccb.jhu.edu/software/stringtie/index.shtml?t=manual

＃https://qiita.com/rouninnomi/items/5441bef2f50780035127

後半文字の大きさがみづらくなってすみません。修正法がわかりませんでした。

Hisat2がerror

mikan03 — Sat, 21 Jul 2018 15:14:13 +0900

Hisat2のindexを作る作業が動きません。このようなエラーがでます。

$ hisat2-build /Mus_musculus/UCSC/mm10/Sequence/WholeGenomeFasta/genome.fa
genome_index
dyld: Symbol not found:
__ZNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEED1Ev
Referenced from: /usr/local/Cellar/hisat2/2.1.0/bin/hisat2-build-s Expected in: /usr/lib/libstdc++.6.0.9.dylib in /usr/local/Cellar/hisat2/2.1.0/bin/hisat2-build-s Abort trap: 6

いろいろ調べて試しましたが、意味がわかりませんでした。問題解決のヒントをご教示いただければ幸いです。 Brew update doctor upgradeは問題なかったです。

configです。

$ brew config
HOMEBREW_VERSION: 1.7.0
ORIGIN: https://github.com/Homebrew/brew
HEAD: 7e8fb9a0f8ab5e841763ec6c7fefa72a8462b594
Last commit: 5 days ago
Core tap ORIGIN: https://github.com/Homebrew/homebrew-core
Core tap HEAD: fa5253faa1d416851bbd3d957f67f5ccbcb15441
Core tap last commit: 6 hours ago
HOMEBREW_PREFIX: /usr/local
CPU: quad-core 64-bit skylake
Homebrew Ruby: 2.3.7 => /usr/local/Homebrew/Library/Homebrew/vendor/portable-ruby/2.3.7/bin/ruby
Clang: 9.1 build 902
Git: 2.18.0 => /usr/local/bin/git
Curl: 7.54.0 => /usr/bin/curl
Java: 10.0.2, 1.8.0_121
macOS: 10.13.6-x86_64
CLT: 9.4.1.0.1.1528165917
Xcode: 9.4.1
XQuartz: 2.7.11 => /opt/X11

どうぞよろしくお願いします。

ChIP-seqデータの比較について

juyoutai4423 — Thu, 15 Mar 2018 17:30:33 +0900

ChIP-seqは定性的であり、定量的な比較はできないと習いました。

しかし、論文の図でヒストンアセチル化等のヒストグラムをよく目にします。

縦軸がtag densityとなっていますが、このtag densityとは何でしょうか？一種の正規化なのでしょうか？

また、この場合tag densityの大小でChIP-seqデータを定量的に比較可能でしょうか？両者を比較してヒストンアセチル化が多い、少ないというのは可能でしょうか？

どなたかご存じの方がいましたらよろしくお願いします。

公共のChIP-seqデータのデータセット間のピーク比較について

takebo — Fri, 12 Jan 2018 18:25:33 +0900

ここ数ヶ月でNGSデータ解析について学んでいます。現在は主に公共データのChIP-Seqデータのピーク検出や比較を行っています。今回新たに利用したいデータセットがあるのですが，その利用について少々悩んでおりますので本サイトに投稿しています。詳しい先生方がいらっしゃいましたらご助言頂けると幸いです。

質問内容は以下の通りです。

異なる機種のIllumina HiSeqを用いて得たSRAデータは互いにデータとして比較することができるのでしょうか。

たとえば，公共ゲノムデータベースGEOに公開されている， GSE86164：Illumina HiSeq 2500を用いて得たChIP-Seqデータ GSE76655：Illumina HiSeq 2000を用いて得たChIP-Seqデータ GSE47043：Illumina HiSeq 2000を用いて得たChIP-SeqデータとRNA-Seqデータは，FASTQ変換，マッピング，ピーク検出などの処理をした後，ピークコールの違いとして互いのデータを比較できるのでしょうか。このように公共ゲノムデータベースのデータセット間での比較について悩んでおります。

どうぞよろしくお願いいたします。

非モデル動物におけるtranscriptome referenceを用いた発現量解析

SOCM — Wed, 10 Jan 2018 19:14:06 +0900

去年DRY解析教本を読んでMacを買い、色々なサイトで勉強させて頂いているバイオインフォ初心者です。この度、かなり行き詰まってしまいましたので、こちらの先生方にご相談させて頂きたく存じます、よろしくお願いいたします。

・やりたいこと

RNA-seqデータ(HiSeq2000 paired-end read)をtranscriptomeのreferenceにマップして発現量解析。最終的にはcummeRubundで各サンプル間の遺伝子発現の違いを表示したいと考えております。

・やったこと

bowtie2でリードをtranscriptome referenceにマッピング、得られたbamファイルをeXpressで定量しFPKMやTPMを算出。しかし各ライブラリーごとの数値を統合するところで手段を失いました。

・問題点

解析対象の動物種おけるGTFファイルが存在しないため（ゲノム情報が整理されていない）、genomeをreferenceとした一連のTophat-cufflinks-cuffdiffのワークフローが使えない。また、transcriptomeをreferenceとして発現量解析をするワークフローに関する情報が乏しい。以上のことから、何が現在もっとも正確なtranscriptomeをreferenceとした発現量解析なのか自身には判断できない状態になっております。

どんな情報でもありがたいです、どなたかアドバイス頂きたく存じます。よろしくお願いいたします。

tophatのエラーについて

SK — Tue, 31 Oct 2017 21:59:50 +0900

dry解析教本という本で自習中の初心者です。大学のMacでは出来たのに家のMacMiniだとエラーが出て、理由が分かりません。

SKMac-mini:expression SK$ tophat -p 1 -G Homo_sapiens/NCBI/build37.2/Annotation/Archives/archive-2014-06-02-13-47-29/Genes/genes.gtf -o tophat_results/ERR266337_P5 Homo_sapiens/NCBI/build37.2/Sequence/Bowtie2Index/genome ERR266337_trim.fastq

[2017-10-31 21:32:10] Beginning TopHat run (v2.1.1)

[2017-10-31 21:32:10] Checking for Bowtie Bowtie version: 2.3.3.1 [2017-10-31 21:32:10] Checking for Bowtie index files (genome).. [2017-10-31 21:32:10] Checking for reference FASTA file [2017-10-31 21:32:10] Generating SAM header for Homo_sapiens/NCBI/build37.2/Sequence/Bowtie2Index/genome Traceback (most recent call last): File "/usr/local/bin/tophat", line 4107, in <module> sys.exit(main()) File "/usr/local/bin/tophat", line 3961, in main params.read_params = check_reads_format(params, reads_list) File "/usr/local/bin/tophat", line 1856, in check_reads_format zf = ZReader(f_name, params) File "/usr/local/bin/tophat", line 1809, in init self.file=open(filename) IOError: [Errno 2] No such file or directory: 'ERR266337_trim.fastq'

ところが今いるexpressionというディレクトリには確かにERR266337_trim.fastqが存在するはずだと思います。

SKMac-mini:expression SK$ ls ERR266335.fastq ERR266347.fastq ERR266351_trim.fastq　 ERR266335_trim.fastq　 ERR266347.fastq.gz FastQC ERR266337.fastq ERR266347_trim.fastq　 Homo_sapiens ERR266337_trim.fastq　 ERR266349.fastq Homo_sapiens_NCBI_build37.2.tar.gz ERR266338.fastq ERR266349_trim.fastq　 tophat_results ERR266338_trim.fastq　 ERR266351.fastq

このエラーはどうして出てきたのでしょう？ご助言がいただければ幸いです。

異なる2系統からF2ゲノムを作成してくれるシミュレーションプログラム

myoshi — Thu, 05 Oct 2017 00:26:46 +0900

実験動物等の異なる2系統からF2ゲノムを作成してくれるシミュレーションプログラムを探しています。両親のfastqファイルから作成するのでも、vcfから作成するのでもいいです。どなたかご存じの方いらっしゃいましたら、ご教授いただけますと助かります。目的はlinkage analysisツールを作ったのですが、そのvalidationです。

heminested PCRの省略

deer — Tue, 03 Oct 2017 18:38:19 +0900

もうひとつ初学者的質問をさせてください。

あるウィルス属をouter primer set(sense A, antisense B)で増幅したのち、その中のウィルス種をinner primer set (sense A, antisense C)で検出する手順が文献にのっています。

属を検出することに興味がなくて種だけ検出したい場合は(sense A, antisense C)でいきなりPCRしてしまってもよいように思いますが、いかがでしょうか？　antisense BとantisenseCは４ntほどかぶっています。

ちなみに1段階目はEx Taq polymeraseの存在下で50μLの反応物を35サイクル(94℃ 30秒　に最初のサイクルだけ2分追加、55℃で30秒、72℃で30秒、最終サイクルに7分追加）。2段階目は、先の2μLにPlaitnum taq polymeraseを加えて、25サイクル（94℃　30秒、最初のサイクルだけ2分追加、55℃で30秒、72℃で30秒、最終サイクルに7分追加）とあります。

2段階を省略して1段階だけにできないものかどうか、ご指導お願いします。

微妙な変異ウィルスの定量方法

deer — Sun, 01 Oct 2017 16:18:14 +0900

お世話になります。実験科学は素人なのでご指導お願いします。（実験は共同研究機関や外注会社にお願いしています）

ウィルス変異株の量（viremia)が病気の重症度と関係があるか調べたいと思っています。

たとえば、ある病原性ウイルスに株１、株２、株３の三つがあるとします。これらの差異は微妙な点突然変異でアミノ酸の違いもわずかであり、それらの違いを区別するPCR primerはまだ設計されていません。おのおのの検体中で株１、株２、株３の「量的」な違いを調べるにはどうすればよろしいでしょうか？ PCRベースの実験ではどうすれば可能でしょうか？

あるいは、株１，２，３に共通のprimerでPCRしておいて（あるいはPCRはなしで）、イルミナにかけたのち、株１，２，３おのおののReferene 配列にmapされるread数を比較してもいいように思います。しかしウィルス量がlog copy/mLのような絶対的な単位ででないのが信頼度の面でどうかなと思いますのと、mapされたread数の多寡がウィルス量の多寡といえるのか自分でもよくわかりません。統計解析するにしてもRead数の%を独立変数としてよいのかどうか。。。また、検体の数が増えてくると処理やコストが（PCRに比して）たいへんになってくるのでは、と危惧しています。

微妙なウィルス変異株の定量方法についてご指導ください。

TMM正規化法について

abcde2017 — Mon, 11 Sep 2017 15:13:55 +0900

次世代シーケンサーからのデータを解析し，発現差異解析を行いたいと思っています。 Aligned BAMという形式のデータをTMM正規化し，発現差の解析をする場合，データに含まれる遺伝子数がTMM正規化により減ってしまうことはありますか？

NCBI GEO ではどのような解析がされているのでしょうか

sunday — Wed, 12 Jul 2017 15:44:09 +0900

NCBI GEO 初学者です。下記につきどなたか詳細をご存知の方がいらっしゃったらご教示願えないでしょうか。

https://goo.gl/wqe5Xd
上記 URL のように NCBI GEO を検索すると遺伝子が 7 個抽出されます (最後のは該当せず)。"disease state"[FINF] により疾患群と対象群とが比較され、違いのあるものが提示されているようです。問題は、この "発現頻度に違いのある遺伝子" の抽出アルゴリズムが良く分かりません。恐らく fold change などを元に出しているのだと考えますが、説明文を見つけることができませんでした。

どうぞ宜しくお願いいたします。

不完全長18Sを用いた系統樹作成に関して

nomat — Tue, 06 Jun 2017 11:50:39 +0900

18Sの配列を用いた系統樹解析を行う場合に、18Sが完全長でない場合にどのように解析するのが一般的なのでしょうか。

現在、自身の持つサンプルの18S配列を用いて近縁種との系統樹を作成しようとしているのですが、自身の持つ18Sが完全長ではなく、短いものだと1.1kb程度となっています。当初は系統樹作成にClustalWを用いたマルチプルアライメント+最尤法での系統樹作成を考えていたのですが、配列長が異なるため、グローバルアライメントでアライメントしてしまうClustalWでは不適切なのでは？と疑問を覚えました。一般的にはこのような場合どのような形で解析を行うのでしょうか。

近縁種の完全長18S配列と特別な処理をせずに解析を実施するのか、それともローカルアライメントにより共通領域を洗い出し、部分領域でのみ比較をして系統樹を作成するのでしょうか。ご教授いただけましたら幸甚と存じます。

cuffmergeの結果について

mlck — Tue, 21 Mar 2017 17:48:07 +0900

御世話になります．初心者なので、DDBJ Sequence Read Archive (DRA) から取得したRNA-seqのヒトサンプルデータ数名分を解析しております．

cufflinksのcuffmergeを実行しました．群内における発現量のばらつきを知りたいのですが、mergedファイル内のgenes.fpkm_tracking（isoforms.fpkm_trackingとどちらが適切かわかりませんが）の結果をみると、FPKM_conf_loとFPKM_conf_hiは同じ値になっています．

これは、うまくmergeできなかったと解釈すればよいでしょうか．

実行コマンドは、cuffmerge -o merged -p 4 -g Homo_sapiens.GRCh38.87.gtf transcripts.gtf.txt　です．

どうぞよろしくお願い致します．

BioMart V0.9を使い倒す　のリンク不具合

KK — Fri, 23 Dec 2016 21:15:42 +0900

発現制御解析・可視化＞遺伝子、タンパク質を機能別に分類したい　＞　BioMart V0.9を使い倒す　のリンクですが、クリックするとその下のDavidの使い方　実践編に飛んでしまいます。(Mac OS Sierra, Firefox 50.1.0) これは使用環境が不適切なためでしょうか、それともリンク切れ（？）なのでしょうか。

よろしくお願いします。

RNA-seqデータを用いたエンリッチメント解析について

lonicera — Mon, 15 Aug 2016 10:54:15 +0900

DAVIDを用いてエンリッチメント解析を行いたいのですが、UploadしたIDが上手く認識されず困っています。

自分が持っているID

RNA-seqにより得られた発現遺伝子のデータに対して、BLASTXを用いて相同性検索を行い、「GI number」と「RefSeq ID」を得ています。

行った動作

「GI number」や「RefSeq ID」を用いてDAVIDによるエンリッチメント解析を行いました。『Select Identifier』の項目ではそれぞれ「PROTEIN_GI_ACCESSION」と「REFSEQ_PROTEIN」を選んでいます。2300のIDを用いて調べた結果、認識されるのは232と非常に少ない結果しか得られませんでした。DAVIDに用意されている「Gene ID conversion」を用い、「AFFYMETRIX」のIDに変換を試みましたが、変換されるIDの数が少なく未だ解析できずにいます。BioMartを用いてIDを変換しようとしましたが、現在一時的にサイトが使えなくなっていました。

みなさんはRNA-seqから得られたデータを用いてエンリッチメント解析を行う際、どのようにしているのでしょうか？また、自分が行った行動の中で不適切な点があれば、ご教授をお願いいたします。

ホモログの検出について

lonicera — Thu, 14 Jul 2016 18:05:18 +0900

RNA-seqを用いて研究を行い始めたものです。右も左もわからない中、ここでの回答は非常助かっております。

現在得られたRNA-seqの発現遺伝子データセットに対して、既存遺伝子のホモログが存在するのかどうかを調べたいと考えています。データセットの量が膨大で、どう調べればよいかが分かりません。

初歩的な質問ではありますが、ご教授して頂けないでしょうか。何卒よろしくお願いいたします。

【追記1】他の生物の遺伝子を対象とし、私の対象とする生物の発達過程でどのような発現をしているかを調査したいと考えています。

回答ありがとうございました。データセットのアノテーション情報を参考にして遺伝子を絞り、調査してきたいと思います。

1000人ゲノムデータのvcfファイルからのDAFの計算方法について

suimye — Tue, 05 Jan 2016 17:17:41 +0900

suimyeです。こんにちは。
表題の件について、ご存知の方是非ご教授いただきたくpostいたしました。
1000人ゲノムデータのVCFファイルを使ってDerived Allele Freq.（以後DAF）を計算したいのですが、調べた結果がどうも納得がいかなくてpostしております。

DAFは、祖先アレル（以後AA）に対して、祖先アレルと一致しないアレルの頻度を計算するものと考えております。

例:

祖先アレル: G

Aさんのアレル: G

Bさんのアレル: A

Cさんのアレル: G

Derived Allele: A

この場合のDAFは、1/3となる。簡単には以上のようなものと思います。

次に、1000人ゲノムのデータのVCFファイルを例にすると、

1 10177 rs367896724 A AC 100 PASS AC=2130;AF=0.425319;AN=5008;N S=2504;DP=103152;EAS_AF=0.3363;AMR_AF=0.3602;AFR_AF=0.4909;EUR_AF=0.4056;SAS_AF=0.49 49;AA=A|||;VT=INDEL

このデータの場合は、AAがAなので、referenceと同じ塩基であり、Derived Alleleとして2塩基の"AC"がコールされています。ヒトの各地域集団ごとのアレル頻度や、ヒト全体でのアレル頻度はAFまたは~_AFで与えられていますが、これはAlternateであるACのアレル頻度なのでDAFはこの頻度をそのまま使えばよいと思っていました。ところが、Biostarなどの幾つかのpostを見ると、

Please read the vcf spec first. AF is the alternate allele frequency. AA is the ancestral allele. If both present, you can use it to get the derived allele frequency by flipping AF.

のようにアレル頻度をフリップして計算しろと書いてあります。このコメントは、そのvariantに祖先アレルがある場合はアレル頻度（AF）の逆、DAF = 1 - AFすれば良いとのことだと思いますが納得できずにいます。
私の拙い考えですと先にも書きましたように、AF値0.425319はreferenceに対するAlternateなので、variantとしてみつかったACのアレル頻度を示しており、DAF=0.425319であるように思います。これは間違いなのでしょうか？くだんのBiostarでは、例にとられたVCFの行には祖先アレルはcallされていませんので、referenceとalternateのどちらが祖先アレルなのか分らず、回答者の「AFをフリップしたらDAF」という回答には釈然としない状況です。もしご存知の方いらっしゃいましたらご教示いただけますと幸甚です。

@suimye

IGVの結果とFPKM値の違いについて

K_Yamashita — Wed, 25 Nov 2015 23:31:59 +0900

NGSデータ解析の初心者です。現在、RNA-seq解析を行っているのですが、Mapping結果をIGVで可視化するとある遺伝子に対してMappingされていることを確認したのですが、その遺伝子の発現量をCuffdiffで計算したところ、FPKM値が0となりました(gene tracking)。リードがMappingされているのにFPKM値がなぜ0となるのかわからず困っています。何か知見があればご教授ください。よろしくお願い致します。

癌のゲノム解析データベースに対する配列検索

Hiroya — Sat, 17 Oct 2015 20:01:21 +0900

NCBIなどに蓄積している癌のゲノムデータベースに対してある特別な配列（人にはない配列です）のインテグレーションを調べたいと考えています。出来ればすべての癌に対して調べたいのですが、難しいなら特定の癌腫に限定しようと思っています。どのような手順を踏むのが最短でできますでしょうか？全てのゲノム配列をダウンロードするのは、かなり大変な作業だと思いますが、ダウンロードする他ないでしょうか？事前の予想では、予めalignmentの取られている配列だと、除外されているのではと考えていますが、現段階では分かりません。

上記の質問は具体的ではなかったので、変更させていただきます。 https://cghub.ucsc.edu/ 上記よりゲノムデータをダウンロードして、Blast等でサーチできる環境を作りたいのですが、具体的にどのような手順を踏むのがベストでしょうか？

GEO からデータのダウンロード：MAX formatとは？

megu — Wed, 02 Sep 2015 16:52:25 +0900

平素お世話になっております。浅学のため記入内容および語彙の不備があるかと思いますが、どうか御容赦頂ければと思います。

●背景：GEOからRNA-seqのデータを取得

NCBI が運営するGEOからRNA-seqのデータを取得したいと考えております。目的とするデータのアクセッション番号はGSE20116です。配列データとしては6ファイルあるようです（例えばGSM515513）。

●質問：MAX formatとは？

サンプル GSM515513を見ますと、以下のような記述がありました。

Data processing AB WT Pipeline. The reads were aligned to the human hg18 (NCBI Build 36) genome. Alignments are provided in MAX format. A MAX specification file is available on the Series GSE20116 record.

上記のMAX formatとは何でしょうか？いろいろ検索してみましたが判然としませんでした。

●困っている現状

実は背景として、DDBJのftpサーバーから上記サンプルに相当するfastqファイルを取得し、tophatにてhg19へのアライメントを試みた経緯があります。ところが、accepted_hits.bamが 1 MB未満となってしまい、明らかに結果がおかしかったです。海外の掲示板も参照し、似たような事例が散見されましたが、解決できず往生しております。GEO のサイトから、マッピング済みのファイルが得られるなら、取得したいと考えています。MAX format なるものが、tophat に適用出来るかはわかりませんが、そもそもMAX formatがなんなのかすら分からない現状です。

些細な点でも構いませんので、どうか御教授頂けないでしょうか？お手数をおかけ致しますが、何卒宜しくお願い致します。

敬具

MacOS X 10.9.5に対するaugustus3.1のインストール時のエラーについて

fatman_2 — Sun, 19 Jul 2015 00:25:37 +0900

お世話になります. 初歩的な質問で申し訳有りませんが, ご教授いただけると幸いです.

ゲノム解読したデータをもとに, Mac OSでab initio遺伝子予測をしようと考えています. そこで, MacOS X 10.9.5のPCに対してaugustus3.1のインストールを試みているのですが, makeがうまくいかず, 以下のようなエラーメッセージが表示されます.

=======================

error: invalid operands to binary expression ('const GBFeature' and 'const GBFeature')

error: invalid operands to binary expression ('const OptionListItem' and 'const OptionListItem')

error: invalid operands to binary expression ('const GBFeature' and 'const GBFeature')

3 warnings and 2 errors generated. make: *** [genbank.o] Error 1

=======================

なにが問題でmakeできないのでしょうか.

よろしくお願いいたします.

RNA seqにおけるRを用いたマッピング時のエラーについて

brown — Mon, 06 Apr 2015 17:58:56 +0900

「Rで塩基配列解析」をみながら、RNAseq解析をしようとしています。Rは初心者で、試しに適当なリファレンス配列とクエリー配列（共にfasta形式）を用いてmappingを試みています。リストファイルを作成後、以下のコマンドによりQuasRを起動してマッピングを行おうとしましたが、以下のようにbowtieにおいて、”引数の長さが 0 です　sh: line 1: 9966 Segmentation fault: 11”というエラーが出てしまい困っております。このような場合の解決策をご存知の方はいらっしゃいますか？

in_f1 <- "mapping1.txt"

in_f2 <- "Ala2.fa"

library(QuasR)

out <- qAlign(in_f1, in_f2)

Creating .fai file for: /Users/ka/Desktop/te/Ala2.fa

alignment files missing - need to: create alignment index for the genome create 1 genomic alignment(s) will start in ..9s..8s..7s..6s..5s..4s..3s..2s..1s

Creating an Rbowtie index for /Users/ka/Desktop/te/Ala2.fa

Finished creating index Testing the compute nodes...OK Loading QuasR on the compute nodes...OK

Available cores: nodeNames ka.biol.sci.u.ac.jp 1

Performing genomic alignments for 1 samples. See progress in the log file: /Users/ka/Desktop/te/QuasR_log_1d53f48759d.txt

以下にエラー checkForRemoteErrors(val) :

one node produced an error: Error on ka.biol.sci.u.ac.jp processing sample /Users/ka/Desktop/te/Ala1.fa : 引数の長さが 0 です sh: line 1: 9966 Segmentation fault: 11

'/Library/Frameworks/R.framework/Versions/3.1/Resources/library/Rbowtie/bowtie'

'/Users/ka/Desktop/te/Ala2.fa.Rbowtie/bowtieIndex' '/Users/ka/Desktop/te/Ala1.fa' -m 1 --best --strata -v 2 -f -S -p 1 '/var/folders/kg/hrs2ndkn029f978fs1_18s140000gn/T//RtmpF7YdOB/Ala1.fa26e25902bc72.sam' 2>&1

下記にsessionInfo() による情報と使用Macのスペックを示します。どうぞよろしくお願い致します。

sessionInfo() R version 3.1.3 (2015-03-09) Platform: x86_64-apple-darwin10.8.0 (64-bit) Running under: OS X 10.8.5 (Mountain Lion)

locale: [1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8

attached base packages: [1] stats4 parallel stats graphics grDevices utils datasets [8] methods base

other attached packages: [1] GenomicAlignments_1.2.2 Rsamtools_1.18.3
[3] Biostrings_2.34.1 XVector_0.6.0
[5] QuasR_1.6.2 Rbowtie_1.6.0
[7] GenomicRanges_1.18.4 GenomeInfoDb_1.2.4
[9] IRanges_2.0.1 S4Vectors_0.4.0
[11] BiocGenerics_0.12.1 BiocInstaller_1.16.2

loaded via a namespace (and not attached): [1] AnnotationDbi_1.28.2 base64enc_0.1-2 BatchJobs_1.6
[4] BBmisc_1.9 Biobase_2.26.0 BiocParallel_1.0.3
[7] biomaRt_2.22.0 bitops_1.0-6 brew_1.0-6
[10] BSgenome_1.34.1 checkmate_1.5.2 codetools_0.2-11
[13] DBI_0.3.1 digest_0.6.8 fail_1.2
[16] foreach_1.4.2 GenomicFeatures_1.18.6 grid_3.1.3
[19] hwriter_1.3.2 iterators_1.0.7 lattice_0.20-31
[22] latticeExtra_0.6-26 RColorBrewer_1.1-2 RCurl_1.95-4.5
[25] RSQLite_1.0.0 rtracklayer_1.26.3 sendmailR_1.2-1
[28] ShortRead_1.24.0 stringr_0.6.2 tools_3.1.3
[31] XML_3.98-1.1 zlibbioc_1.12.0

Mac Proスペックプロセッサ 2 x 2.66 GHz 6-Core Intel Xeon メモリ 32 GB 1333 MHz DDR3 ECC ソフトウェア OS X 10.8.5（12F45）

RNA-seqのAT bias

junya — Mon, 06 Apr 2015 17:38:15 +0900

MiSeqを使ってmouse macrophage系細胞のRNA-seqを始めたのですが、出てきたデータの塩基組成がA,T＝32%、G,C＝18%と偏ります。原因に心当たりのある方、アドバイス頂けないでしょうか？

実験系の概略は以下のとおりです。

Total RNA 1 ugを精製
TruSeq Stranded mRNA Sample Prep Kitでlibrary構築
Sequencing Kit V3 (150cycle)でpaired-end readsを取得

よろしくお願いいたします。

MISOアノテーションについて

Kent_allow — Tue, 17 Feb 2015 17:41:46 +0900

質問させていただきます。

現在、選択的スプラシング解析ツール"MISO"を用いた解析を行なっております。
これについて数点問題があり、どなたかご教示いただければと思った次第です。

1)
現在使用しているアノテーションおよびインデックスは、UCSC genome browserよりDLしたgtfファイルを、rnaseqlibに含まれる"gff make annotation.py"というスクリプトを用いて手前でビルドしています。
この中身を確認した所、過去にASが報告されているisoformに関するエキソンが入っておりませんでした。
UCSC genome browserよりDLしたgtfには、これらのisoformが含まれていることを確認しています。
この様な事例をご存知の方がいましたら、原因・対策について教えていただければ幸いです。

2)
上記に関する可能性の一つとして、一部のisoformのevent typeによるものが考えられます。
報告分子のevent typeがALEである可能性があるのですが、gff make annotation.pyによりビルドされたアノテーションにはALE のものは含まれません。
そこで、(MISO Exon-centric annotation ver. 1のように)ALE、AFEを含めたアノテーションをビルドする方法をご存知の方がいましたらご教示いただければ幸いです。

以上ご存知の方がいましたら、どうぞよろしくお願い致します。

NCBIのMicrobial Nucleotide BLASTをLocal環境に構築したいのですが、

TY — Fri, 13 Feb 2015 01:20:14 +0900

NCBI BLASTのトップページに行くと、一番上にBLAST Assembled Genomesという項目があると思います。 Human, Mouseなど代表的な生物種が並んで、最後のMicrobesをクリックすると、Microbial Nucleotide BLAST という検索フォームが開きます。

Choose Search Setの項目の中のDatabaseとして、「Representative genomes only」と「All genomes」が選べ、「All genomes」を選ぶとさらに

Complete genomes
Draft genomes
Complete plasmid
Complete bacteriophages

が選べるようになっています。これらのデータベースをLocal環境にダウンロードしたいのですが、 FTPサイト (ftp://ftp.ncbi.nlm.nih.gov/blast/db/) を見ても、どれがそのデータにあたるのかが分からず、困っています。

更新日時と、結果画面で表示される「Included: Archaea (taxid:2157), Bacteria(taxid:2)」という表現から推測して、「refseq_genomic.[数字].tar.gz」というデータベースをtaxidで絞り込んで使用しているようなのですが、これをどうBLASTコマンドで再現すればいいかが分かりません。

統合TVの「Local BLAST の使い方 2011」シリーズ（全2回）

第1回導入・準備編
第2回検索実行・オプション編

は拝見したのですが、taxidでの絞り込みはなかったので質問させて頂きました。どうぞよろしくお願い致します。

RNA-seq で得られるFPKM のグラフの作り方

ara — Mon, 26 Jan 2015 16:50:51 +0900

いつも勉強させて頂いております。使用している語彙が不適切かもしれませんがどうかご容赦頂きたく思います。

お忙しいところ畏れ入りますがご教授頂けますと幸いです。

★背景

細胞A と細胞B について、Biological triplicates (n=3)　を用意しRNA-sequencing を実施しました。つまり、細胞Aについて3回分のRNA-seqのデータ、細胞Bについて3回分のRNA-seqのデータがあります。現在、TopHat、Cufflinks, Cuffcompare, Cuffdiff で得られたファイルが手元にあります。

★質問

これらのデータを使って、細胞Aと細胞Bに発現する遺伝子（一種類のみ）のFPKMを比較したグラフを書きたいと考えております。恥ずかしながら浅学のため、以下の点についてご意見を頂けないでしょうか？

① グラフには、Cuffdiff を実行した際に得られる「genes.fpkm_tracking」中のFPKMを使用して宜しいでしょうか？

② ①のFPKM でグラフを作成した場合、いわゆるエラーバーには、どのファイルのどの値を用いれば良いのでしょうか？

③ それとも各RNA-seq データずつFPKM を算出し、それをエクセル上で平均値、標準偏差を求めてグラフを作成するべきでしょうか？

質問が初歩的で大変申し訳ありません。Cuffdiff は「total reads 数を考量し triplicates をまとめたFPKM」が得られる、と理解しておりますが、そのFPKM でグラフを作るのが適切なことなのかどうか、そしてエラーバーのための値（標準偏差）は？と往生しております。

もちろん、そもそも根本的に間違っていることも覚悟しております。その場合、恥ずかしながらいわゆるスタンダードなグラフ作成方法をご教授頂けますと助かります。

お手数をおかけ致しますどうか何卒宜しくお願い致します。

未回答質問

アミノ酸配列でサイトごとに計算される保存度は何を意味しているのか？

参考書に載っていた内容

マルチプルアラインメントされたアミノ酸配列で、サイトごとに保存度を「以下の式」で計算する意義について知りたいです。

保存度を計算する理由について、私の考察

考察に対する反論

よって、私の考察とは異なる目的があるのだろう。という結論にたどり着きました。ご回答いただけると幸いです。

pyBedtoolsによるベン図の作成

BAMファイルのソートについて

ngsplotdbのインストールについて

DESeq2でtranscripts IDとGene IDを紐付けする方法

追記：Rでの作業を書きます。

公式ページに従って入力。transcritp_count_matrix.csvはstringtieで生成

合計6未満(N=3で２郡間の比較）の微弱発現を除去

PHENO_DATA.txtでサンプル割り当て。#PHENO_DATA.txtの中身

Genotype

1 WT

2 WT

3 WT

4 KO

5 KO

6 KO

確認1

[1] TRUE がでるのでOK

[1] TRUE がでるのでOK

DEseq2に流し込み。

Gene IDをつけるためrtracklayarうごかす。

Stringtieでmergeしたファイルをいれる。

stringtieの場合はtypeがexonかtranscript のようなので、exonにする。

紐付けされたリストを取り出しておく。

CSVへの書き出しで終わるはず。

参考ページ(どうもありがとうございます。）

https://ncrna.jp/blog/item/388-deseq2-ggplot2

http://ccb.jhu.edu/software/stringtie/index.shtml?t=manual

後半文字の大きさがみづらくなってすみません。修正法がわかりませんでした。

Hisat2がerror

ChIP-seqデータの比較について

公共のChIP-seqデータのデータセット間のピーク比較について

非モデル動物におけるtranscriptome referenceを用いた発現量解析

tophatのエラーについて

[2017-10-31 21:32:10] Beginning TopHat run (v2.1.1)

異なる2系統からF2ゲノムを作成してくれるシミュレーションプログラム

heminested PCRの省略

微妙な変異ウィルスの定量方法

TMM正規化法について

NCBI GEO ではどのような解析がされているのでしょうか

不完全長18Sを用いた系統樹作成に関して

cuffmergeの結果について

BioMart V0.9を使い倒す のリンク不具合

RNA-seqデータを用いたエンリッチメント解析について

ホモログの検出について

1000人ゲノムデータのvcfファイルからのDAFの計算方法について

IGVの結果とFPKM値の違いについて

癌のゲノム解析データベースに対する配列検索

GEO からデータのダウンロード：MAX formatとは？

MacOS X 10.9.5に対するaugustus3.1のインストール時のエラーについて

RNA seqにおけるRを用いたマッピング時のエラーについて

RNA-seqのAT bias

MISOアノテーションについて

NCBIのMicrobial Nucleotide BLASTをLocal環境に構築したいのですが、

RNA-seq で得られるFPKM のグラフの作り方

BioMart V0.9を使い倒す　のリンク不具合