すべての質問

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

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

考察に対する反論

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

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

アミノ酸配列の保存度算出に伴う重み付けスコアの意義について

study_bioinfo — Mon, 16 Sep 2019 01:57:26 +0900

はじめまして、初めて使用させていただきます。よろしくお願いします。

1,疑問点

・以下の論文でValder氏の提唱している式(37式)を学びましたが、その式中でHenikoff氏の考案した重み付けがなぜ使用されているのかわからない。 (参照:https://www.ncbi.nlm.nih.gov/pubmed/12112692)

疑問を持った理由

→例えば、100本のマルチプルアラインメントされた相同配列のうち、99本が90%程度の一致度を示し、残りの1本が20%程度の一致度を示しているような場合を考える。重み付けの計算式によれば、20%の一致度を示す配列の情報を保存度に大きな影響として与えようとしていることがわかる。なぜ、このような計算方法をとる必要があるのか？99本の保存されていると考えられるサイトを以って、全ての配列は90%程度保存されていると判断してはいけないのはなぜか？と考えたため。

---論文を読み込めない方のために、使用している式について以下に画像をキャプチャしました---

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

2,重み付けの計算方法(https://gyazo.com/7907a5beddc2a1ee6348e310f9a5d78b)

※iは配列の数を表し、Nは配列の総数を表している。xはポジションのことを表している。

→お分かりになる方がいるだろうという希望的観測のもと質問をいたしますので、計算方法の詳細については記述いたしません。必要な場合は追記いたしますので、おしらせください。

3,スコア行列から取得したアミノ酸ペアの値に対する計算(https://gyazo.com/60242cf27db66929a5490dc0ed9b1558)

4,スコア行列=Blosum62(https://gyazo.com/7c8c0bb28db10ba26d07bdf6e40aed57)

2,私の理解度について

Valder氏の式については、"加重平均"を計算していることは理解しております。配列ごとに重み付けをした値を分母と分子にとり、分子についてはBlosum62のスコア行列から得られた値を配列同士の重み付けの積と、さらに積をとって計算しています。分子、分母共に、各項が1を超えることは絶対にないように計算されているので、保存度は0~1の間で算出される。

3,疑問点解決のために模索したこと

・Henikoff氏の重み付けを使用する理由は、「taxonomic bias」を逃れるために導入していることを以下の参考文献から学びましたが、このバイアスについて具体例を元に詳しく説明されている文献等々を見つけることができず、具体的に記している文献を知りたいです。 (参考文献:タンパク質機能解析のためのバイオインフォマティクス P.67)

4,求めている回答

・taxonomic biasを理解することが、重み付けを式に組み込む理由として正しいのか？

→妥当である場合、taxonomic biasが発生してしまうような場合を、具体例を含んで解説している書籍、論文をご紹介いただきたいです。

→妥当でない場合、何を理解することで重み付けをする理由を説明することができるかをご指南いただきたいです。

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

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

Tophatとcufflinksがインストールできない

mikan03 — Thu, 19 Jul 2018 15:29:31 +0900

こんにちは

Tophatとcufflinksがインストールできません。レベルは初心者なので、わけもわからずやっている状態ですが、 DRY解析教本を読みながらコツコツやっています。いろいろやってみましたが、どうにもなりませんでした。

$brew install tophat Error: No available formula with the name "tophat" ==> Searching for a previously deleted formula (in the last month)... Error: No previously deleted formula found. ==> Searching for similarly named formulae... Error: No similarly named formulae found. ==> Searching taps... ==> Searching taps on GitHub... Error: No formulae found in taps.

$ brew install cufflinks Error: No available formula with the name "cufflinks" ==> Searching for a previously deleted formula (in the last month)... Error: No previously deleted formula found. ==> Searching for similarly named formulae... Error: No similarly named formulae found. ==> Searching taps... ==> Searching taps on GitHub... Error: No formulae found in taps.

$brew tap brewsci/science ==> Tapping brewsci/science Cloning into '/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science'... remote: Counting objects: 468, done. remote: Compressing objects: 100% (465/465), done. remote: Total 468 (delta 1), reused 153 (delta 1), pack-reused 0 Receiving objects: 100% (468/468), 429.98 KiB | 305.00 KiB/s, done. Resolving deltas: 100% (1/1), done. Error: Invalid formula: /usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/ome-common.rb uninitialized constant #<class:0x00000001032840c0>::MinimumMacOSRequirement Did you mean? MaximumMacOSRequirement Error: Invalid formula: /usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/ome-files.rb uninitialized constant #<class:0x00000001023892a0>::MinimumMacOSRequirement Did you mean? MaximumMacOSRequirement Error: Invalid formula: /usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/ome-xml.rb uninitialized constant #<class:0x00000001022e3a08>::MinimumMacOSRequirement Did you mean? MaximumMacOSRequirement Error: Invalid formula: /usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/osgearth.rb uninitialized constant #<class:0x0000000103086fe8>::MinimumMacOSRequirement Did you mean? MaximumMacOSRequirement Error: Cannot tap brewsci/science: invalid syntax in tap!

$ brew config
HOMEBREW_VERSION: 1.7.0
ORIGIN: https://github.com/Homebrew/brew
HEAD: 7e8fb9a0f8ab5e841763ec6c7fefa72a8462b594
Last commit: 3 days ago
Core tap ORIGIN: https://github.com/Homebrew/homebrew-core
Core tap HEAD: ea3fec61c268ccbbe4172a2f1cbfd08760025cf0
Core tap last commit: 51 minutes 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

$ brew list
bam     eigen       ilmbase     numpy       samtools
bedtools    fastqc      jpeg        opencv      sqlite
blast       ffmpeg      lame        openexr     stringtie
boost       gdbm        libidn2     openssl     tbb
bowtie      gettext     libpng      python      wget
bowtie2     git     libtiff     python@2    x264
bwa     hisat2      libunistring    readline    xvid
curl        htslib      lmdb        rmtrash     xz

brewのupdate, upgrade, doctorはしました。問題なかったです。問題を解決するにあたって何を質問したら良いかもわかっていない状態です。

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

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を含めたアノテーションをビルドする方法をご存知の方がいましたらご教示いただければ幸いです。

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

すべての質問

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

参考書に載っていた内容

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

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

考察に対する反論

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

アミノ酸配列の保存度算出に伴う重み付けスコアの意義について

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

Tophatとcufflinksがインストールできない

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アノテーションについて

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