十二 24

TCGA数据下载大全

Posted on 2015年12月24日 by ulwvfje

并不是所有的数据都能下载，很多数据需要有权限才能下载的！！！

首先我们可以根据TCGA的文章来下载数据：

https://tcga-data.nci.nih.gov/docs/publications/

总共也就几十篇文章，都是发表在大杂志上面的。

每篇文章都会提供数据的打包下载，例如：

The molecular taxonomy of primary prostate cancer
Cell Volume 163 Issue 4: p1011-1025 Read the full article
Portal Publication Site and Associated Data Files

Comprehensive Molecular Characterization of Papillary Renal Cell Carcinoma
NEJM. Published on line on Nov 4th, 2015 Read the full article
Portal Publication Site and Associated Data Files

那个portal链接点击进去，就可以看到所有的下载链接了！

https://tcga-data.nci.nih.gov/docs/publications/prad_2015/

这是根据文章来分类打包好的数据！

同时也可以通过其它数据接口来下载

Tools for Exploring Data and Analyses

Broad Institute FireBrowse portal, The Broad Institute
cBioPortal for Cancer Genomics, Memorial Sloan-Kettering Cancer Center
TCGA Batch Effects, MD Anderson Cancer Center
Regulome Explorer, Institute for Systems Biology
Next-Generation Clustered Heat Maps, MD Anderson Cancer Center

TCGA Data Portal

这几个接口都挺好用的：

每个接口都有使用教程：http://firebrowse.org/tutorial/FireBrowse-Tutorial.pdf

非常详细！！！而且还专门写了软件接口：https://confluence.broadinstitute.org/display/GDAC/Download

或者写了R的接口：http://www.cbioportal.org/cgds_r.jsp

一般都推荐用TCGA自己的数据接口：https://tcga-data.nci.nih.gov/tcga/

里面对所有的样本都进行了统计

通过https://tcga-data.nci.nih.gov/tcga/dataAccessMatrix.htm可以进行定制化的数据下载！

里面有很多TCGA自定义的名词：

Data Levels and Data Types

https://tcga-data.nci.nih.gov/docs/dictionary/ 可以看到所有名词的解释：

数据的种类如下：

还记得以前看到一篇TCGA自己的关于胃癌的文章，发表在nature上面，文章涉及到了TCGA的各个方面的分析，所以附件PDF竟然有133页！！！

文章地址是：http://www.nature.com/nature/journal/v513/n7517/full/nature13480.html

pdf附件地址是：http://www.nature.com/nature/journal/v513/n7517/extref/nature13480-s1.pdf

所有的数据地址：https://tcga-data.nci.nih.gov/docs/publications/stad_2014/

设计到的295个sample的list在：https://tcga-data.nci.nih.gov/tcgafiles/ftp_auth/distro_ftpusers/anonymous/other/publications/stad_2014/STAD.Sample_Barcodes.txt

包含的其它数据有：

十二 23

看看Y染色体上面的基因在测序数据里的覆盖度和测序深度

Posted on 2015年12月23日 by ulwvfje

作用：可以检测别人是否把自己的样本搞混，也可以看看测序是否分布均匀！

首先，我们要拿到Y染色体上面的基因的坐标信息！

因为我们的是hg19，所以我们要下载hg19的基因信息！

wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/refGene.txt.gz

我们首先解析refGene文件，找到chrY的unique基因！

这四列分别是：chromosome/start/end/gene_symbol

4程序如下：

[perl]
open FH,"/home/jmzeng/hg19/chrY.gene.special.position" or die "file error !!!";
while(<FH>){
 chomp;
 @F=split;
 foreach ($F[1]..$F[2]){
 $h{$_}=$F[3];
 }
 $length{$F[3]}=$F[2]-$F[1]+1;
}
close FH;
open FH,$ARGV[0];
while(<FH>){
 chomp;
 @F=split;
 next unless $F[0] eq 'chrY';
 next if $F[2]<20;
 if (exists $h{$F[1]}){
 $count{$h{$F[1]}}++ ;
 }else{
 $count{'other'}++ ;
 }
}
close FH;
print "$_\t$length{$_}\t$count{$_}\n" foreach sort keys %count;</pre>
</div>
<div>[/perl]

对一个男性样本，结果会如下：

gene/length/pos

AMELY	8111	1269
BCORP1	47724	689
CSPG4P1Y	3799	538
DAZ1	69739	762
DAZ2	71901	228
DAZ3	73222	233
DAZ4	73222	540
DDX3Y	12825	3654
EIF1AY	17445	929
FAM224A	4295	82
FAM224B	4293	85
GOLGA2P3Y	4866	68
GYG2P1	15476	547
HSFY2	42277	3950
KDM5D	39526	7425
NLGN4Y	319396	3872
PCDH11Y	105374	6627
PRKY	107577	1390
PRORY	3388	735
RBMY1B	14451	232
RBMY1D	14411	117
RBMY1E	14410	157
RBMY1J	14407	65
RBMY2EP	6416	27
RBMY2FP	7348	419
RPS4Y1	25376	1856
RPS4Y2	24966	1831
SRY	888	703
TBL1Y	180999	3231
TGIF2LY	958	808
TMSB4Y	2457	534
TSPY4	132211	1525
TTTY14	205048	394
TTTY4C	36811	39
TTTY9A	9317	580
TXLNGY	23067	1968
USP9Y	159610	10508
UTY	232293	6670
VCY	742	291
XKRY2	1582	980
ZFY	47437	3125
other		100328

对女性样本，结果会如下；

NLGN4Y	319396	575
PCDH11Y	105374	1643
PRKY	107577	82
TGIF2LY	958	191
TTTY14	205048	139
other		54297

从结果可以看出来，很多基因都是y染色体特有的，这个结果是表明我们的测序非常棒

十二 22

根据chrY独有区域的覆盖度及测序深度来判断性别

Posted on 2015年12月22日 by ulwvfje

这个也是基于bam文件来的，判断chrY独有区域的覆盖度及测序深度

首先下载chrY独有区域的记录文件，https://www.familytreedna.com/documents/bigy_targets.txt

然后用samtools depth来统计测序深度，samtools depth $i |grep 'chr[XY]'

depth统计结果文件如下：

mzeng@ubuntu:/home/jmzeng/gender_determination$ head Sample3.depth
chrX   60085   1
chrX   60086   1
chrX   60087   1
chrX   60088   1
chrX   60089   1
chrX   60090   1
chrX   60091   1
chrX   60092   1
chrX   60093   1
chrX   60094   1

然后我随便写了一个脚本来对测序深度文件进行再统计，统计覆盖度及测序深度

[perl]
open FH,"bigy_targets.txt";
while(<FH>){
chomp;
@F=split;
$all+=$F[2]-$F[1]+1;
foreach ($F[1]..$F[2]){
$h{$_}=1;
}
}
close FH;
open FH,$ARGV[0];
while(<FH>){
chomp;
@F=split;
next unless $F[0] eq 'chrY';
if (exists $h{$F[1]}){
$pos++;
$depth+=$F[2];
}
}
close FH;
$average=$depth/$pos;
$coverage=$pos/$all;
print "$pos\t$average\t$coverage\n" ;

[/perl]

这样对那三个样本结果如下：

可以看到只有sample4，是覆盖率极低的，而且记录到的pos位点也特别少，所以她是女性！

这里测序深度没有意义。

chrY独有区域文献：https://www.familytreedna.com/learn/wp-content/uploads/2014/08/BIG_Y_WhitePager.pdf

区域记录于：https://www.familytreedna.com/documents/bigy_targets.txt

另外的文献：http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2748900/#bib8

一个很有趣的博客：http://davetang.org/muse/2013/09/07/creating-a-coverage-plot-in-r/

try the GATK "DepthOfCoverage" ?http://www.broadinstitute.org/gsa/wiki/index.php/Depth_of_Coverage_v3.0

or you can run 'samtools pileup' and calculate the mean value for the coverage column.

十二 22

根据X,Y染色体比对上的reads数来判断性别

Posted on 2015年12月22日 by ulwvfje

针对高通量测序数据，包括WGS，WES

我这里主要讲根据bam文件里面的chrX和chrY的reas比例来判断性别，大家可以自己处理数据得到bam文件！

主要是读取bam文件，选择chrX上面的记录，统计genotype即可：

也可以统计测序深度，samtools depth $i |grep 'chr[XY]'

如果chrX和chrY的reads比例超过一定值，比如50倍，就判定为女性！

脚本很简单，就统计bam文件的第三列就可以啦，第三列就是染色体信息

samtools view $i | perl -alne '{$h{$F[2]}++}END{print "$_\t$h{$_}" foreach sort keys %h}' >$out.chromosome.stat

对于Sample3，

chrX 1233061

chrY 140506

对于Sample4，很明显，这个是女性！

chrX 2734815

chrY 51860

对于Sample5

chrX 1329302

chrY 156663

十二 22

根据X染色体的snp的纯和率来判断性别

Posted on 2015年12月22日 by ulwvfje

针对高通量测序数据，包括WGS，WES，甚至snp6芯片也行。

我这里主要讲根据vcf文件里面的chrX的纯合率来判断性别，大家可以自己处理数据得到vcf文件！

主要是读取vcf文件，选择chrX上面的记录，统计genotype即可：

我这里拿之前的自闭症项目数据来举例子：

根据数据提供者的信息，3-4-5分别就是孩子、父亲、母亲，统计chrX的snp的的纯合和杂合的比例，代码很简单

vcf文件一般第一列是染色体，第6列是质量，第10列是基因型已经测序深度相关信息

cat Sample5.gatk.UG.vcf |perl -alne '{next unless $F[0] eq "chrX";next unless $F[5]>30;$h{(split/:/,$F[9])[0]}++}END{print "$_\t$h{$_}" foreach keys %h}'

如果纯合的snp是杂合的倍数超过一定阈值，比如4倍，就可以判断是男性。

对于Sample3来说，很明显，是男孩，因为X染色体都是纯合突变

0/1 391
1/1 2463
2/2 1
1/2 32
0/2 1
对于Sample4来说，很明显，应该是母亲，证明之前别人给我的信息有误
1/1 3559
1/2 27
0/1 1835
0/2 5
那么Sample5很明显就是父亲咯
1/1 2626
0/1 356
1/2 22

十二 21

用重抽样+主成分方法来做富集分析

Posted on 2015年12月21日 by ulwvfje

之前我们用超几何分布检验的方法做了富集分析，使用的是GSE63067.diffexp.NASH-normal.txt的logFC的绝对值大于0.5，并且P-value小雨0.05的基因作为差异基因来检验kegg的pathway的富集情况

结果是这样的

我们接下来用另外一种方法来做富集分析，顺便检验一下，是不是超几何分布统计检验的富集分析方法就是最好的呢？

这种方法是-重抽样+主成分分析

大概的原理是，比如对上图中的，04380这条pathway来说，总共有128个基因，那么我从原来的表达矩阵里面随机抽取128个基因的表达矩阵做主成分分析，并且抽取一千次，每次主成分分析都可以得到第一主成分的贡献度值。那么，当我并不是随机抽取的时候，我就抽04380这条pathway的128个基因，也做主成分分析，并且计算得到第一主成分分析的重要性值。我们看看这个值，跟随机抽1000次得到的值差别大不大。

这时候就需要用到表达矩阵啦！

setwd("D:\\my_tutorial\\补\\用limma包对芯片数据做差异分析")

exprSet=read.table("GSE63067_series_matrix.txt.gz",comment.char = "!",stringsAsFactors=F,header=T)

rownames(exprSet)=exprSet[,1]

exprSet=exprSet[,-1]

我们根据ncbi里面对GSE63067的介绍可以知道，对应NASH和normal的样本的ID号，就可以提取我们需要的表达矩阵

把前面两属于Steatosis的样本去掉即可，exprSet=exprSet[,-c(1:2)]

然后再把芯片探针的id转换成entrez id

exprSet=exprSet[,-c(1:2)]

library(hgu133plus2.db)

library(annotate)

platformDB="hgu133plus2.db";

probeset <- rownames(exprSet)

rowMeans <- rowMeans(exprSet)

EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))

match_row=aggregate(rowMeans,by=list(EGID),max)

colnames(match_row)=c("EGID","rowMeans")

dat=data.frame(EGID,rowMeans,probeset)

tmp_prob=merge(dat,match_row,by=c("EGID","rowMeans"))

relevantProbesets=as.character(tmp_prob$probeset)

length(relevantProbesets) #hgu133plus2.db 20156

exprSet=exprSet[relevantProbesets,]

EGID_name=as.numeric(lookUp(relevantProbesets, platformDB, "ENTREZID"))

rownames(exprSet)=as.character(EGID_name)

d=exprSet

最后得到表达矩阵表格

我们首先得到1000次随机挑选128个基因的表达矩阵的主成分分析，第一主成分贡献度值。

gene128=sapply(1:1000,function(y) {

dat=t(d[sample(row.names(d), 128, replace=TRUE), ]);

round(100*summary(fast.prcomp(dat))$importance[2,1],2)

}

)

很快就能得到结果，可以看到数据如下

> summary(gene128)

Min. 1st Qu. Median Mean 3rd Qu. Max.

19.1 25.8 28.8 29.8 32.5 59.7

那么接下来我们挑选这个04380这条pathway特有128个基因来算第一主成分贡献度值

path_04380_gene=intersect(rownames(d),as.character(Path2GeneID[['04380']]))

dat=t(d[path_04380_gene,]);

round(100*summary(fast.prcomp(dat))$importance[2,1],2)

得到的值是38.83，然后看看我们的这个38.83在之前随机得到的1000个数里面是否正常，就按照正态分布检验来算

1-pnorm((38.83-mean(gene128))/sd(gene128))

[1] 0.0625

可以看到已经非常显著的不正常了，可以说明这条通路被富集了。

至少说明超几何分布检验的方法得到的富集分析结果跟我们这次的重抽样+主成分分析结果是一致的，当然，也有不一致的，不然就不用发明一种新的方法了。

如果写一个循环同样可以检验所有的通路，但是这样就不需要事先准备好差异基因啦！！！这是这个分析方法的特点！

十二 21

主成分分析略讲

Posted on 2015年12月21日 by ulwvfje

主成分分析是为了简化变量的个数。
我这里不涉及到任何高级统计知识来简单讲解一下主成分分析，首先我们用下面的代码随机创造一个矩阵：

options(digits = 2)
x=c(rnorm(5),rnorm(5)+4)
y=3*c(rnorm(5),rnorm(5)+4)
dat=rbind(x,y,a=0.1*x,b=0.2*x,c=0.3*x,o=0.1*y,p=0.2*y,q=0.3*y)
colnames(dat)=paste('s',1:10,sep="")
dat
library(gmodels)
pca=fast.prcomp(t(dat))
pca
summary(pca)$importance
biplot(pca, cex=c(1.3, 1.2));

Continue reading →

十二 15

用超几何分布检验做富集分析

Posted on 2015年12月15日 by ulwvfje

我们可以直接使用R的bioconductor里面的一个包，GOstats里面的函数来做超几何分布检验，看看每条pathway是否会富集

我们直接读取用limma包做好的差异分析结果

setwd("D:\\my_tutorial\\补\\用limma包对芯片数据做差异分析")

DEG=read.table("GSE63067.diffexp.NASH-normal.txt",stringsAsFactors = F)

View(DEG)

我们挑选logFC的绝对值大于0.5，并且P-value小雨0.05的基因作为差异基因，并且转换成entrezID

probeset=rownames(DEG[abs(DEG[,1])>0.5 & DEG[,4]<0.05,])

library(hgu133plus2.db)

library(annotate)

platformDB="hgu133plus2.db";

EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))

length(unique(EGID))

#[1] 775

diff_gene_list <- unique(EGID)

这样我们的到来775个差异基因的一个list

首先我们直接使用R的bioconductor里面的一个包，GOstats里面的函数来做超几何分布检验，看看每条pathway是否会富集

library(GOstats)

library(org.Hs.eg.db)

#then do kegg pathway enrichment !

hyperG.params = new("KEGGHyperGParams", geneIds=diff_gene_list, universeGeneIds=NULL, annotation="org.Hs.eg.db",

categoryName="KEGG", pvalueCutoff=1, testDirection = "over")

KEGG.hyperG.results = hyperGTest(hyperG.params);

htmlReport(KEGG.hyperG.results, file="kegg.enrichment.html", summary.args=list("htmlLinks"=TRUE))

结果如下：

但是这样我们就忽略了其中的原理，我们不知道这些数据是如何算出来的，只是由别人写好的包得到了结果罢了。

事实上，这个包的这个hyperGTest函数无法就是包装了一个超几何分布检验而已。

如果我们了解了其中的统计学原理，我们完全可以写成一个自建的函数来实现同样的功能。

超几何分布很简单，球分成黑白两色，数量已知，那么你随机抽有限个球，应该抽多少白球的问题！

公式就是 exp_count=n*M/N

然后你实际上抽了多少白球，就可以计算一个概率值！

换算成通路的富集概念就是，总共有多少基因，你的通路有多少基因，你的通路被抽中了多少基因（在差异基因里面属于你的通路的基因），这样的数据就足够算出上面表格里面所有的数据啦！

tmp=toTable(org.Hs.egPATH)

GeneID2Path=tapply(tmp[,2],as.factor(tmp[,1]),function(x) x)

Path2GeneID=tapply(tmp[,1],as.factor(tmp[,2]),function(x) x)

#phyper(k-1,M, N-M, n, lower.tail=F)

#n*M/N

diff_gene_has_path=intersect(diff_gene_list,names(GeneID2Path))

n=length(diff_gene_has_path) #321 # 这里算出你总共抽取了多少个球

N=length(GeneID2Path) #5870 ##这里算出你总共有多少个球（这里是错的，有多少个球取决于背景基因！一般是两万个）

options(digits = 4)

for (i in names(Path2GeneID)){

M=length(Path2GeneID[[i]]) ##这个算出你的所有的球里面，白球有多少个

exp_count=n*M/N ###这里算出你抽取的球里面应该多多少个是白色

k=0 ##这个k是你实际上抽取了多少个白球

for (j in diff_gene_has_path){

if (i %in% GeneID2Path[[j]]) k=k+1

}

OddsRatio=k/exp_count

p=phyper(k-1,M, N-M, n, lower.tail=F) ##根据你实际上抽取的白球个数，就能算出富集概率啦！

print(paste(i,p,OddsRatio,exp_count,k,M,sep=" "))

}

随便检查一下，就知道结果是一模一样的！

十二 15

下载所有的酶的信息，并且解析好

Posted on 2015年12月15日 by ulwvfje

之前我提到过kegg数据库里面的有些pathway下面没有对应任何基因，当时我还在奇怪，怎么会有这样的通路呢？

然后，我随机挑选了其中一条通路（hsa01000），进行查看，发现正好是所有的酶的信息。

好奇怪，不明白为什么kegg要列出所有的酶信息

http://www.genome.jp/kegg-bin/get_htext?hsa01000

下载htext格式的酶的信息

798K Dec 15 2015 hsa01000.keg

查看文件，发现也是层级非常清楚的结构，D已经是最底级别的酶了，而E对应的基因是属于该酶的。

简单统计一下，发现跟酶相关的基因有3688个，而最底级别的酶有5811个，应该会持续更新的

如果想做成kegg那样的基因与酶对应表格，也是非常简单的！

国际系统分类法按酶促反应类型,将酶分成六个大类:

1、氧化还原酶类(oxidoreductases) 催化底物进行氧化还原反应的酶类。包括电子或氢的转移以及分子氧参加的反应。常见的有脱氢酶、氧化酶、还原酶和过氧化物酶等。

2、转移酶类(transferases) 催化底物进行某些基团转移或交换的酶类,如甲基转移酶、氨基转移酶、转硫酶等。

3、水解酶类(hydrolases)催化底物进行水解反应的酶类。如淀粉酶、粮糖苷酶、蛋白酶等。

4、裂解酶类(lyases)或裂合酶类(synthases) 催化底物通过非水解途径移去一个基团形成双键或其逆反应的酶类,如脱水酶、脱羧酸酶、醛缩酶等。如果催化底物进行逆反应,使其中一底物失去双键,两底物间形成新的化学键,此时为裂合酶类。

5、异构酶类(isomerases)催化各种同分异构体、几何异构体或光学异构体间相互转换的酶类。如异构酶、消旋酶等。

6、连接酶类（ligases）或合成酶类(synthetases)催化两分子底物连接成一个分子化合物的酶类。

上述六大类酶用EC(enzyme commission)加1.2.3.4.5.6编号表示,再按酶所催化的化学键和参加反应的基团,将酶大类再进一步分成亚类和亚-亚类,最后为该酶在这亚-亚类中的排序。如α淀粉酶的国际系统分类纩号为:EC3.2.1.1

EC3——Hydrolases 水解酶类

EC3.2——Glycosylases 转葡糖基酶亚类

EC3.2.l——Glycosidases 糖苷酶亚亚类i.e.enzymes hydmlyzing O-and S-glycosyl compound即能水解O-和S-糖基化合物

EC3.2.1.1 Alpha-amylase, α-淀粉酶

值得注意的是,即使是同一名称和EC编号,但来自不同的物种或不同的组织和细胞的同一种酸,如来自动物胰脏、麦芽等和枯草杆菌BF7658的α-淀粉酶等,它们的一级结构或反应机制可解不同,它们虽然都能催化淀扮的水解反应,但有不同的活力和最适合反应条件。

可以按照酶在国际分类编号或其推荐名,从酶手册(Enzyme Handbook)、酶数据库中检索到该酶的结构、特性、活力测定和Km值等有用信息。著名的手册和数据库有:

手册:

1、Schomburg,M.Salzmann and D.Stephan:Enzyme Handbook 10 Volumes

2、美国Worthington Biochemical Corporation:Enzyme Manual

(http://www.worthington-biochem.com/index/manual.htm/)

数据库:

l、德国BRENDA：Enzyme Database(http://www.brenda wnzymes.org)

2、Swissprot:EXPASYENZYME Enzyme nomenclature database (http://www.expasy.org/enzyme/)

3、IntEnz:Integrated relational Enzyme database (http://www.ebi.ac.uk/mtenz)

十二 13

可怕的中国教育（转）

Posted on 2015年12月13日 by ulwvfje

无意中看到这篇文章，感觉挺有趣的！

我不知道国外的教育是咋样的，但是文中描述的中国教育的试题都是千真万确的，我也经历过。

我没看过有钱人家的小孩，不知道那些花了几百万买个几平米的学区房的那些人的小孩是不是也要经受这样的教育摧残？

我还有个问题，我是如何出淤泥而不染的呢？

随便瞎扯，大家自己看这篇文章吧！

可怕的中国教育，聪明伶俐进去呆若木鸡出来

作者：葛小琴

　　来源：杂文月刊（文摘版）

　　侄子在读高二，考了一道历史题：成吉思汗的继承人窝阔台，公元哪一年死?最远打到哪里?答不出来，我帮他查找资料，所以到现在我都记得，是打到现在的匈牙利附近。

　　在一次偶然的机会，我发现美国世界史这道题目不是这样考的。它的题目是这样的：成吉思汗的继承人窝阔台，当初如果没有死，欧洲会发生什么变化?试从经济、政治、社会三方面分析。

　　有个学生是这样回答的：这位蒙古领导人如果当初没有死，那么可怕的黑死病，就不会被带到欧洲去，后来才知道那个东西是老鼠身上的跳蚤引起的鼠疫。但是六百多年前，黑死病在欧洲猖獗的时候，谁晓得这个叫做鼠疫?如果没有黑死病，神父跟修女就不会死亡。神父跟修女如果没有死亡，就不会怀疑上帝的存在。如果没有怀疑上帝的存在，就不会有意大利弗罗伦斯的文艺复兴?

　　如果没有文艺复兴，西班牙、南欧就不会强大，西班牙无敌舰队就不可能建立。如果西班牙、意大利不够强大，盎格鲁—撒克逊会提早200年强大，日耳曼会控制中欧，奥匈帝国就不可能存在。

　　教师一看“棒，分析得好。”但他们没有分数，只有等级A。其实这种题目老师是没有标准答案的，可是大家都要思考。

　　不久前，我去了趟日本，日本总是和我们在历史问题上产生纠葛，所以我在日本很注意高中生的教科书。

　　他们的教师给高中生布置了这样一道题：日本跟中国100年打一次仗，19世纪打了日清战争(即甲午战争)，20世纪打了一场日中战争(即抗日战争)，21世纪如果日本跟中国开火，你认为大概是什么时候?可能的远因和近因在哪里?如果日本赢了，是赢在什么地方?输了是输在什么条件上?分析之。

　　其中有个高中生是这样分析的：我们跟中国很可能在台湾回到中国以后，有一场激战。台湾如果回到中国，中国会把基隆与高雄封锁，台湾海峡就会变成中国的内海，我们的油轮就统统走右边，走基隆和高雄的右边。这样，会增加日本的运油成本。我们的石油从波斯湾出来跨过印度洋，穿过马六甲海峡，上中国南海，跨台湾海峡进东海到日本海，这是石油生命线，中国政府如果把台湾海峡封锁起来，我们的货轮一定要从那里经过，我们的主力舰和驱逐舰就会出动，中国海军一看到日本出兵，马上就会上场，就开打。

　　按照判断，公元2015年至2020年之间，这场战争可能爆发。所以，我们现在就要做对华抗战的准备。

　　我看其他学生的判断，也都是中国跟日本的摩擦会从东海从台湾海峡开始，时间判断是2015年至2020年之间。

　　这种题目和答案都太可怕了。

　　撇开政治因素来看这道题，我们的历史教育就很有问题。翻开我们的教科书，题目是这样出的：甲午战争是哪一年爆发的？签订的叫什么条约？割让多少土地？赔偿多少银两？每个学生都努力做答案。结果我们一天到晚研究什么时候割让辽东半岛，什么时候丢了台湾、澎湖、赔偿二万银两，1894年爆发甲午战争、 1895年签订马关条约，背得滚瓜烂熟，都是一大堆枯燥无味的数字。

　　那又怎么样，反正都赔了嘛，银两都给了嘛，最主要的是将来可能会怎样。

　　人家是在培养能力，而我们是在灌输知识，这是值得深思的部分！

　　看外面的教育，再看我们的教育？

　　老妈去参加我侄子的家长会，带回了两套侄子的考试试卷，我很好奇，拿过来看了现在小学生的试卷后，我震惊了！这是什么狗屁教育？这样的教育有希望吗？下面给大家详细说说我看到了什么。

　　侄子在本市某著名小学读书，有这么几道题。

　　一个春天的夜晚，一个久别家乡的人，望着皎洁的月光不禁思念起了故乡，于是吟起了一首诗：( )，( )。

　　我看到侄子答的是：举头望明月，低头思故乡。但后面是一把大大的×，我就奇怪了，我也是想到的这两句。好奇地问侄子，这个不对??那答案是什么?侄子说标准答案是：春风又绿江南岸，明月何时照我还。哎，这就奇怪了，因为是个春天的夜晚，就要是这句有春风的？？？要这个思念故乡的人不是江南的，是不可能说出春风又绿江南岸这句话的!!!举头望明月，低头思故乡应该更准确。再扯远点，思念故乡，一千个人可以吟一千句不一样的诗，这个也可以有标准答案的么?

　　接下来是默写，题目是：我们学过《桂林山水》一文，请将下面句子默写下来。然后就是整段的要默写，这有什么用？死记硬背别人的文字有什么用?

　　还有个题目，《匆匆》这篇课文，是现代着名作家朱自清先生写的，同学们都很喜欢这篇散文，你能把自己最喜欢，印象最深刻的一句写下来吗？我侄子写的是：我的日子滴在时间的流里，没有声音，也没有影子。后面一把好大的×。标准答案竟然是：但是，聪明的你告诉我，我们的日子为什么一去不复返呢？这就更奇怪了，一篇文章，你可以喜欢这句，我可以喜欢那句，难道最喜欢的一句话也要统一么?为什么“我的日子滴在时间的流里，没有声音，也没有影子。”这句不能喜欢?就一定要喜欢“但是，聪明的你告诉我，我们的日子为什么一去不复返呢?”这句？我觉得这个题目应该是“你能把老师最喜欢，印象最深刻的一句写下来吗?”才对!

　　再看别的试卷，更莫名其妙了，比如请说出阿拉伯数字的来历，是哪个国家创造的？侄子不知道，问我，我也不知道。我只好去搜一下，才知道是古印度人发明的。莫非我吃块猪肉，还一定得知道它是哪个养猪场养出来的?

　　最后有个题目让我彻底崩溃了：请用一句话说明π的含义。侄子回答π的含义是圆周率。竟然打的是×，这就奇怪了，正好我老婆大学说读的是理科，我马上问她，π是什么意思，她说圆周率啊。两个人狂汗，问了侄子半天，标准答案大概是，π是一个在数学及物理学领域普遍存在的数学常数……

　　如果你也觉得这种教育很无耻，就请转发吧，让更多的人来参与呼吁改变，为了孩子为了国家的未来 …… 别让孩子聪明伶俐地进去呆若木鸡地出来！

十二 12

下载最新版GO，并且解析好

Posted on 2015年12月12日 by ulwvfje

首先要明白，需要下载什么？

要下载四万多条GO记录的详细信息（http://purl.obolibrary.org/obo/go/go-basic.obo）

要下载GO与GO之间的关系（http://archive.geneontology.org/latest-termdb/go_daily-termdb-tables.tar.gz）

要下载GO与基因之间的对应关系！（物种）（ftp://ftp.ncbi.nlm.nih.gov/gene/DATA）

去官网！

http://geneontology.org/page/download-ontology

grep '\[Term\]' go-basic.obo |wc

43992 43992 307944

版本的区别!刚才我们下载的GO共有43992条，而以前的版本才38804条

GO与GO之间的关系

对应关系也在更新

> as.list(GOBPPARENTS['GO:0000042'])

$`GO:0000042`

is_a is_a is_a is_a

"GO:0000301" "GO:0006605" "GO:0016482" "GO:0072600"

library(org.Hs.eg.db)

library(GO.db)

> tmp=toTable(org.Hs.egGO) ##这个只包括基因与最直接的go的对应关系

> dim(tmp)

[1] 213101 4

> tmp2=toTable(org.Hs.egGO2ALLEGS) #这个是所有的基因与go的对应关系

> dim(tmp2)

[1] 2218968 4

基因与GO的对应关系也在更新

grep '^9606' gene2go |wc -l ### ##这个只包括基因与最直接的go的对应关系

269063

ftp://ftp.informatics.jax.org/pub/reports/index.html#go

http://www.ebi.ac.uk/QuickGO

ftp://ftp.ncbi.nlm.nih.gov/gene/DATA

ftp://ftp.informatics.jax.org/pub/reports/index.html#go

十二 11

关于limma包差异分析结果的logFC解释

Posted on 2015年12月11日 by ulwvfje

首先，我们要明白，limma接受的输入参数就是一个表达矩阵，而且是log后的表达矩阵（以2为底）。

那么最后计算得到的logFC这一列的值，其实就是输入的表达矩阵中case一组的平均表达量减去control一组的平均表达量的值，那么就会有正负之分，代表了case相当于control组来说，该基因是上调还是下调。

我之前总是有疑问，明明是case一组的平均表达量和control一组的平均表达量差值呀，跟log foldchange没有什么关系呀。

后来，我终于想通了，因为我们输入的是log后的表达矩阵，那么case一组的平均表达量和control一组的平均表达量都是log了的，那么它们的差值其实就是log的foldchange

首先，我们要理解foldchange的意义，如果case是平均表达量是8，control是2，那么foldchange就是4，logFC就是2咯

那么在limma包里面，输入的时候case的平均表达量被log后是3，control是1，那么差值是2，就是说logFC就是2。

这不是巧合，只是一个很简单的数学公式log(x/y)=log(x)-log(y)

十二 11

用excel表格做差异分析

Posted on 2015年12月11日 by ulwvfje

其实主要要讲的不是用excel来做差异分析，只是想讲清楚差异分析的原理，用excel可视化的操作可能会更方便理解，而且想告诉大家，其实生物信息学分析，本来就很简单的，那么多软件，只有你理解了原理，你自己就能写出来的！

首先，还是得到表达矩阵，下面绿色的样本是NASH组，蓝色的样本是normal组

我们进行差异分析，很简单，就是看两组的表达值，是否差异，而检验的方法就是T检验。

=AVERAGE(D2:L2) ##求NASH组的平均表达量

=AVERAGE(M2:S2) ###求normal的平均表达量

=T2-U2 ##计算得到logFOLDchange值

=AVERAGE(D2:S2) ###得到所有样本的平均表达量

=T.TEST(D2:L2,M2:T2,2,3) ###用T检验得到两个组的表达量的差异显著程度。

简单检查几个值就可以看到跟limma包得到的结果差不多。

十二 11

用limma包对芯片数据做差异分析

Posted on 2015年12月11日 by ulwvfje

下载该R语言包，然后看说明书，需要自己做好三个数据（表达矩阵，分组矩阵，差异比较矩阵），总共三个步骤（lmFit,eBayes,topTable）就可以啦

首先做第一个数据，基因表达矩阵！

自己在NCBI里面可以查到下载地址，然后用R语言读取即可

exprSet=read.table("GSE63067_series_matrix.txt.gz",comment.char = "!",stringsAsFactors=F,header=T)

rownames(exprSet)=exprSet[,1]

exprSet=exprSet[,-1]

然后做好分组矩阵，如下

然后做好，差异比较矩阵，就是说明你想把那些组拿起来做差异分析，如下

最后输出结果：

我进行了6次比较，所以会输出6次比较结果

最后打开差异结果，解读，说明书如下！

忒

在我的github有完整代码

十二 08

下载最新版的KEGG信息，并且解析好

Posted on 2015年12月8日 by ulwvfje

打开官网：http://www.genome.jp/kegg-bin/get_htext?hsa00001+3101

http://www.genome.jp/kegg-bin/get_htext#A1 （这个好像打不开）

可以在里面找到下载链接

下载得到文本文件，可以看到里面的结构层次非常清楚，

C开头的就是kegg的pathway的ID所在行，D开头的就是属于它的kegg的所有的基因

A,B是kegg的分类，总共是6个大类，42个小类

grep ^A hsa00001.keg

AMetabolism

AGenetic Information Processing

AEnvironmental Information Processing

ACellular Processes

AOrganismal Systems

AHuman Diseases

也可以看到，到目前为止（2015年12月8日20:26:57），共有343个kegg的pathway信息啦

接下来我们就把这个信息解析一下：

perl -alne '{if(/^C/){/PATH:hsa(\d+)/;$kegg=$1}else{print "$kegg\t$F[1]" if /^D/ and $kegg;}}' hsa00001.keg >kegg2gene.txt

这样就得到了

但是我发现了一个问题，有些通路竟然是没有基因的，我不是很明白为什么？

C 04030 G protein-coupled receptors [BR:hsa04030]

C 01020 Enzyme-linked receptors [BR:hsa01020]

C 04050 Cytokine receptors [BR:hsa04050]

C 03310 Nuclear receptors [BR:hsa03310]

C 04040 Ion channels [BR:hsa04040]

C 04031 GTP-binding proteins [BR:hsa04031]

那我们来看看kegg数据库更新的情况吧。

首先我们看org.Hs.eg.db这个R包里面自带的数据

Date for KEGG data: 2011-Mar15

org.Hs.egPATH has 5869 entrez genes and 229 pathways

2015年八月我用的时候是 6901 entrez genes and 295 pathways

现在是299个通路，6992个基因

所以这个更新其实很缓慢的，所以大家还在用DAVID这种网络工具做kegg的富集分析结果也差不大！

详细信息见http://www.genome.jp/kegg/pathway.html

更新信息见：http://www.genome.jp/kegg/docs/upd_map.html

十二 01

我的生物信息学视频上线啦

Posted on 2015年12月1日 by ulwvfje

虽然优酷比较坑，但是他的受众多，大家可以去http://i.youku.com/trainee 上面找到我的全部视频，免费观看！

视频列表是http://www.bio-info-trainee.com/tmp/tutorial/video_list.html ，能顺利点击的链接代表视频录制完毕。

这次视频课程的大纲是http://www.bio-info-trainee.com/tmp/tutorial/syllabus.htm，但是很有可能会修改！

声明
大家好，欢迎观看由生信菜鸟团举办的生物信息学公益视频课程，我是主讲人Jimmy。
本课程仅面向那些生物出身却想转生物信息学数据分析的同学，其它人均不在考虑范围。
本课程在中国大陆录制，会尽量遵循中国大陆法律，如有法律问题，请联系我的律师，谢谢！
我会尽量保证视频中知识点的准确无误，如有任何讲解不当之处，欢迎批评指正！
发邮件联系我（jmzeng1314@outlook.com）即可，没什么必要不要QQ或者微信找我。

还有，本人不差钱，所以视频均免费的，录制完毕后我会放出百度云共享，现在还在测试阶段。

请不要无缘无故的批评我，我不会服务任何人，所以不可能有人是我的顾客，我不需要对你好，爱看不看！

十一 15

hapmap计划资源收集

Posted on 2015年11月15日 by ulwvfje

官网是：http://hapmap.ncbi.nlm.nih.gov/index.html.en

里面有教程：http://hapmap.ncbi.nlm.nih.gov/tutorials.html.en

所有的数据都放在ncbi上面：ftp://ftp.ncbi.nlm.nih.gov/hapmap/

现在一般用这个计划的数据主要是拿自己得到的突变数据来跟这个hapmap计划的人种突变数据对比。

有芯片数据，也有WES和WGS数据，随着时间的推进，平台也在更新：

Jul 07 2009 00:00    Directory affy100k
Mar 05 2010 00:00    Directory affy500k
Jun 02 2010 00:00    Directory hapmap3_affy6.0

当然，数据也在更新

Jul 07 2009 00:00    Directory 2005-03_phaseI
Dec 03 2009 00:00    Directory 2005-11_phaseII
Jul 07 2009 00:00    Directory 2007-03
Jul 07 2009 00:00    Directory 2008-03
Jul 07 2009 00:00    Directory 2008-07_phaseIII
Jul 07 2009 00:00    Directory 2008-10_phaseII
Jul 07 2009 00:00    Directory 2009-01_phaseIII
Jul 07 2009 00:00    Directory 2009-02_phaseII+III
Aug 18 2010 00:00    Directory 2010-05_phaseIII
Sep 19 2010 00:00    Directory 2010-08_phaseII+III

数据都被整合好了：

Bulk data
- Genotypes: Individual genotype data submitted to the DCC to date. Phase 3 data is available in PLINK format and HapMap format.
- Frequencies: Allele & genotype frequencies compiled from genotyping data submitted to the DCC to date. These have also been submitted to dbSNP and should be available in the next dbSNP build.
- LD Data: Linkage disequilibrium properties D', LOD , R² compiled from the genotype data to date
- Phasing Data: Phasing data generated using the PHASE software, compiled from the genotype data to date.
- Allocated SNPs: dbSNP reference SNP clusters that have been picked and prioritized for genotyping according to several criteria (see info on how SNPs were selected). The file 00README contains per-chromosome SNP counts and further details.
- CNV Genotypes: CNV data from HapMap3 samples.
- Recombination rates and Hotspots: Recombination rates and hotspots compiled from the genotyping data.
- SNP assays: Details about assays submitted to the DCC to date. PCR primers, extension probes etc., specific to each genotyping platform.
- Perlegen amplicons: Details for mapping Perlegen amplicons to HapMap assayLSID. For primer sequences, see Perlegen's Long Range PCR Amplicon data.
- Raw data: Raw signal intensity data from HapMap genotypes. Currently includes data from Affymetrix GeneChip 100k and 500k Mapping Arrays.
- Inferred genotypes: Genotypes inferred using the method of Burdick et al. Nat Genet 38:1002-4.
- Mitochondrial and chrY haplogroups: Classification of phase I HapMap samples into mtDNA and chrY haplogroups. The distribution shown in Table 4 of the HapMap phase I paper (Nat Genet 38:1002-4) corresponds to unrelated parents in each one of the populations analyzed.

同时也发了很多篇文章：

The International HapMap Consortium. Integrating common and rare genetic variation in diverse human populations.
Nature 467, 52-58. 2010. [Abstract] [PDF] [Supplementary information]
The International HapMap Consortium. A second generation human haplotype map of over 3.1 million SNPs.
Nature 449, 851-861. 2007. [Abstract] [PDF] [Supplementary information]
The International HapMap Consortium. A Haplotype Map of the Human Genome.
Nature 437, 1299-1320. 2005. [Abstract] [PDF] [Supplementary information]
The International HapMap Consortium. The International HapMap Project.
Nature 426, 789-796. 2003. [Abstract] [PDF] [Supplementary information]
The International HapMap Consortium. Integrating Ethics and Science in the International HapMap Project.
Nature Reviews Genetics 5, 467 -475. 2004. [Abstract] [PDF]
Thorisson, G.A., Smith, A.V., Krishnan, L., and Stein, L.D. The International HapMap Project Web site.
Genome Research,15:1591-1593. 2005. [Abstract] [PDF]

十一 05

点突变详解

Posted on 2015年11月5日 by ulwvfje

DNA分子中某一个碱基为另一种碱基置换，导致DNA碱基序列异常，是基因突变的一种类型。可分为转换和颠换两类。转换(transitions)是同类碱基的置换（AT→GC及GC→AT）,颠换(transversions) 是不同类碱基的置换（AT→TA或CG,GC→CG或TA。

DNA substitution mutations are of two types. Transitions are interchanges of two-ring purines (A G) or of one-ring pyrimidines (C T): they therefore involve bases of similar shape. Transversions are interchanges of purine for pyrimidine bases, which therefore involve exchange of one-ring and two-ring structures.

我们在分析driver mutation的时候会区分各种点突变:

1. CpG transitions
2. CpG transversions
3. C:G transitions
4. C:G transversions
5. A:T transitions
6. A:T transversions

那么，我们有64种密码子，每种密码子都会有9种突变可能，我们如何得到一个所有的突变可能的分类并且打分表格呢？

类似于下面这样的表格：共576行！！！

head category.acgt
AAA>AAT 2 A T 6
AAA>AAC 2 A C 6
AAA>AAG 2 A G 5
AAA>ATA 1 A T 6
AAA>ACA 1 A C 6
AAA>AGA 1 A G 5
AAA>TAA 0 A T 6
AAA>CAA 0 A C 6
AAA>GAA 0 A G 5
AAT>AAA 2 T A 6

tail category.acgt
GGC>GGG 2 C G 2
GGG>AGG 0 G A 3
GGG>TGG 0 G T 4
GGG>CGG 0 G C 4
GGG>GAG 1 G A 3
GGG>GTG 1 G T 4
GGG>GCG 1 G C 4
GGG>GGA 2 G A 3
GGG>GGT 2 G T 4
GGG>GGC 2 G C 4

我本来以为这是一件很简单的事情，写起来，才发现好麻烦

里面用到的一个函数如下：就是判断突变属于上述六种的哪一种！

参考：https://www.mun.ca/biology/scarr/Transitions_vs_Transversions.html

https://en.wikipedia.org/wiki/Transversion

http://www.uvm.edu/~cgep/Education/Mutations.html

点突变，也称作单碱基替换（single base substitution），指由单个碱基改变发生的突变。

可以分为转换（transitions）和颠换（transversions）两类。

转换：嘌呤和嘌呤之间的替换，或嘧啶和嘧啶之间的替换。

颠换：嘌呤和嘧啶之间的替换。

方便理解下面再附上一张示意图，如下：

十一 05

使用mutsig软件来找驱动基因

Posted on 2015年11月5日 by ulwvfje

从数以万计的突变里面找到driver mutation这个课题很大，里面的软件我接触的就有十几个了，但是我尝试了其中几个，总是无法运行成功，不知道为什么，终于今天成功了一个，就是mutsig软件！其实关于突变数据找driver mutation ，台湾一个大学做了一个数据库DriverDB http://ngs.ym.edu.tw/driverdb/：还因此发了一篇文章：http://nar.oxfordjournals.org/content/early/2013/11/07/nar.gkt1025.full.pdf，挺不错的！

关于driver mutation的理论最近也进化了很多，算是比较完善了吧，但是我一直没时间静下心来好好补充理论知识，很多软件，都只是用过，很多数据，也只是处理了一下，不知道为什么要去做，╮(╯▽╰)╭扯远了，开始谈这个软件吧！

mutsig软件是broadinstitute出品的，所以可靠性非常好咯，来源于一篇nature文章：http://www.nature.com/nature/journal/v505/n7484/full/nature12912.html，而该软件的地址是：http://www.broadinstitute.org/cancer/cga/mutsig_run 需要简单注册才能下载的。

该nature文章是这样描述这个软件的优点的：We used the most recent version of the MutSig suite of tools, which looks for three independent signals: highmutational burden relative to background expectation, accounting for heterogeneity; clustering of mutations within the gene; and enrichment of mutations in evolutionarily conserved sites. Wecombined the significance levels (P values) fromeach test to obtain a single significance level per gene (Methods).

这个软件需要安装matlab环境才能使用，所以我前面就写了教程，如何安装！http://www.bio-info-trainee.com/?p=1166

如果已经安装好了matlab环境，那么直接下载这个软件就可以使用了，软件解压就OK拉，而且人家还提供了测试文件！

软件下载后，解压可以看到里面的一个脚本，软件说明书写的非常简单，当然，使用这个软件也的确非常简单：

run_MutSigCV.sh <path_to_MCR> mutations.maf coverage.txt covariates.txt output.txt 即可，其中所有的数据都是可以下载的，

运行完了测试数据，就证明你的软件安装没有问题啦！如果你只有突变数据的maf格式，maf格式可以参考：https://www.biostars.org/p/69222/ ，也可以使用该软件：如下

run_MutSigCV.sh <path_to_MCR> my_mutations.maf exome_full192.coverage.txt gene.covariates.txt my_results mutation_type_dictionary_file.txt chr_files_hg19

上面三个zip文件，都是可以在mutsig软件官网找到下载链接的，是必须下载的！使用很简单，就一个命令即可，但是把你的vcf突变数据做成该软件需要的maf格式，是一个难题！

十一 05

在linux系统里面安装matlab运行环境mcr

Posted on 2015年11月5日 by ulwvfje

matlab毕竟是收费软件，而且是有界面的。所以搞生物信息的都用R和linux替代了，但是很多高大上的单位，比如大名鼎鼎的broadinstitute，是用matlab的，所以他们开发的程序也会以matlab代码的形式发布。但是考虑到大多研究者用不起matlab，或者不会用，所以就用linux系统里面安装matlab运行环境来解决这个问题，我们仍然可以把人家写的matlab程序，在linux命令行下面，当做一个脚本来运行！

比如，这次我就需要用broadinstitute的一个软件：Mutsig，找cancer driver gene的，http://www.broadinstitute.org/cancer/cga/mutsig_run，但是我看了说明才发现，它是用matlab写的，所以我要想在我的服务器用，就必须按照安装matlab运行环境，在官网可以下载：http://www.mathworks.com/products/compiler/mcr/

我这里选择的是R2013a (8.1)，下载之后解压是这样的，压缩包约四百多M

然后直接在解压后的目录里面运行那个install即可，然后如果你的linux可以传送图像，那么就会想安装windows软件一样方便！如果你的linux是纯粹的命令行，那么，就需要一步步的命令行交互，选择安装地址，等等来安装了。

记住你安装之后，会显示一些环境变量给你，请千万要记住，然后自己去修改自己的环境变量，如果你忘记了，就需要搜索来解决环境变量的问题啦！安装之后是这样的：

请记住你的安装目录，以后你运行其它matlab相关的程序，都需要把这个安装目录，当做一个参数传给你的其它程序的！！！

如果你没有设置环境变量，就会出各种各样的错误，用下面这个脚本可以设置

其中MCRROOT一般是$path/biosoft/matlab_running/v81/ 这样的东西，请务必注意，LD_LIBRARY_PATH非常重要，非常重要，非常重要！！！！

MCRROOT=$1
echo ---
LD_LIBRARY_PATH=.:${MCRROOT}/runtime/glnxa64 ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRROOT}/bin/glnxa64 ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRROOT}/sys/os/glnxa64;
MCRJRE=${MCRROOT}/sys/java/jre/glnxa64/jre/lib/amd64 ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/native_threads ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/server ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/client ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE} ;
XAPPLRESDIR=${MCRROOT}/X11/app-defaults ;
export LD_LIBRARY_PATH;
export XAPPLRESDIR;
echo LD_LIBRARY_PATH is ${LD_LIBRARY_PATH};

如果你没有设置正确，那么会报一下的错误！

error while loading shared libraries: libmwmclmcrrt.so.8.1: cannot o pen shared object file: No such file or directory

error while loading shared libraries: libmwlaunchermain.so: cannot o pen shared object file: No such file or directory

等等！！！

Page 98 of 109« First ‹ Previous 95 96 979899 100 101 Next ›Last »

2024年7月
一	二	三	四	五	六	日
« 四
1	2	3	4	5	6	7
8	9	10	11	12	13	14
15	16	17	18	19	20	21
22	23	24	25	26	27	28
29	30	31