01

脚本作业-解读NCBI的ftp里面关于人的一些基因信息

为了感谢大家对我博客的关注,我在这里发布一个作业,适合菜鸟做的。里面有十几个类似的问题,大家可以下载数据自行处理,如果是问这些问题,我优先回答!

NCBI的ftp里面关于人的一些基因信息

我在NCBI的ftp服务器里面下载了这些数据,时间是2015年,大多是hg19系列的,文件名如下:

CDS.fa 这个是ensembl中人的CDS碱基序列文件,hg38

entrez2go.gene 这个是有go注释的基因情况,有一万八的基因都有go注释

entrez2name.gene 这个是NCBI的entrez ID号对应着基因名的文件

entrez2pubmed.gene 这个是NCBI的entrez ID号对应着该基因发表过的文章的ID号

entrez2refseq2ensembl.gene 这个是NCBI的entrez ID号对应着基因名的refseq的ID号和ensembl数据库的ID号

human_gene_info这个是基因的详细信息,包括基因的起始终止点坐标等等

Protein.fa 这个是ensembl中人的蛋白的氨基酸序列文件,有十万多个蛋白hg38

ref2ensembl.txt  这个是基因名的refseq的ID号和ensembl数据库的ID号

自行去NCBI的ftp服务器里面下载这些数据。

然后好好熟悉这些数据信息,回答一下几个问题:

人总的基因有多少个,它们分别分布在哪些染色体上面,基因的转录本分布情况如何,基因的长度分布如何,基因的外显子个数如何。

CD分子的基因有多少个,它们分别分布在哪些染色体上面,基因的转录本分布情况如何,基因的长度分布如何,基因的外显子个数如何。它们有没有氨基酸偏好性??

MHC系列基因信息?CCL系列基因信息如何?CXCL系列信息如何?或者你感兴趣的基因家族信息?

现在研究最热门的基因是什么?发表文章最多的前十个基因是什么?

基因长度情况如何?最长的基因多长?最短的基因多少bp,可靠吗?

蛋白质长度情况如何?

每条染色体的基因分别情况?基因在染色体那个地方分别最多?

请用图形展示你的结论!!!

 

如果你能回答以上问题,证明你的脚本水平不错了。

如果找不到我,看旁边的公告,加入生信菜鸟群,我就在里面!!!

01

CHIP-seq第三讲之使用MACS软件寻找peaks

在使用Bowtie比对于完Chip-Seq的结果后,就需要用到MACS或者ERANGE来找出峰所在的位置了。但是由于ERANGE的设置比较复杂,所以最为流行的还是MACS。

一.首先安装MACS软件

MACS有两个版本,分别是MACS14和MACS2。MACS2在很多方面都对MACS14做了重大改进,但目前还在测试阶段。我们依然以MACS14为例进行说明。

MACS软件的下载地址在wget https://codeload.github.com/taoliu/MACS/zip/master

这是一个python软件,有152M,已经算是很大了!所以需要按照安装python的方法来安装它!但是,好像这个是最新版的,我们还是用1.4版本吧

wget http://github.com/downloads/taoliu/MACS/MACS-1.4.2-1.tar.gz

其实它的readme已经把这个软件的各种安装使用方法讲的很清楚了。

https://github.com/taoliu/MACS/blob/master/README.rst

MACS软件的具体原理,大家去看文献,或者参考这篇文章

http://www.plob.org/2014/05/08/7227.html

很简单的一个python命令即可安装该软件python setup.py install --user

CHIP-seq第三讲之使用MACS软件寻找peaks752

二.然后准备该软件所需要的数据

是我们在前两篇文章中提到的数据

CHIP-seq第三讲之使用MACS软件寻找peaks786

三.接着运行MACS的命令

/home/jmzeng/.local/bin/macs14 -t Xu_WT_rep2_BAF155.fastq.trimmed.single.bam \

> -c Xu_WT_rep2_Input.fastq.trimmed.single.bam \

> -f BAM -g hs --bw 300 -w -S -n Xu_WT_rep2

CHIP-seq第三讲之使用MACS软件寻找peaks974

四.最后解读一下结果

CHIP-seq第三讲之使用MACS软件寻找peaks987

56K Apr 30 21:54 Xu_WT_rep2_model.r

5.5K Apr 30 22:21 Xu_WT_rep2_negative_peaks.xls

783K Apr 30 22:21 Xu_WT_rep2_peaks.bed

865K Apr 30 22:21 Xu_WT_rep2_peaks.xls

766K Apr 30 22:21 Xu_WT_rep2_summits.bed

唉,反正这也不是我的课题,懒得解释这些结果啦,等后来有机会再慢慢玩吧

 

 

参考 http://www.plob.org/2014/05/08/7227.html

附录:我们现在来了解如何设置参数。

参考自 http://www.plob.org/2014/01/26/7118.html

-t TFILE, –treatment=TFILE 输入文件名

-c CFILE, –control=CFILE 输入阴对文件名

-n NAME, –name=NAME 输入出文件名前缀

-f FORMAT, –format=FORMAT 输入文件格式,默认值为AUTO,可选的值为”BEG”,”ELAND”,”ELANDMULTI”,”ELANDMULTIPER”,”ELANDEXPORT”,”SAM”,”BAM”,”BOWTIE”等。

-g GSIZE, –gsize=GSIZE 比对模板大小。格式可以是:1.0e+9,或者1000000000,也可以缩写:’hs’ for 人类 (2.7e9), ‘mm’ for 大鼠(1.87e9), ‘ce’ for 线虫 (9e7) and ‘dm’ for 果蝇 (1.2e8), 默认值:hs

-s TSIZE, –tsize=TSIZE 设置为短序列的长度,默认值为25

-p PVALE, –pvalue=PVALUE 非峰可能性截取值,默认值为1e-5,这个值不能大太,超过0.9的话,可能无法输出正确的结果

-m MFOLD, –mfold=MFOLD 峰值高度相对于本底的比值,默认值为10,30。也就是说,最低值不能少于10,但比值超过30也不认为它是正常的一个峰。一般而言,低值设置为10是一个很好的区分点。如果这个值还是无法得到满意的结果,那么可以设置得更低,但最好还是使用–nomodel参数,使–nomodel设置为True,然后再传递–shiftsize及–bw参数给MACS。–shiftsize默认值为100,而–bw的默认值为300。

–diag 生成完整报表,会包括是否为真峰的可能性,但会严重拖累运算速度。

 

30

人的CD分子基因信息简介!

CD分子吧,它是Clusters of Differentiation的简写,是指一组分化抗原的家族,目前该家族已经有CD1——CD350甚至更多的成员.他们分布于T细胞等免疫细胞表面,参与免疫细胞各种表达,其中有整合素、受体、配体等蛋白分子,在免疫应答反应中参与识别、粘附和信号转导等功能.

我这里简单讲讲如何整理它们的基因信息,首先从NCBI里面下载的人的gene_info文件,然后通过脚本来查找CD分子信息。

perl -alne '{if (/\tCD\d+/ or /CD\d+\|/ ) {print}}' human_gene_info >CD.info

cut -f 2-5 CD.info >CD.table

再根据CD分子的排序把我们的信息重新排序

perl -alne '{/CD(\d+\w)/;$hash{$1}=$_}END{print $hash{$_} foreach sort {$a <=> $b}keys %hash}' CD.table >CD.table.sort

然后我发现了一个很有趣的问题,它们都是负义链上面的基因!

 

entrez ID gene symbol 正负链
911 CD1C - BDCA1|CD1|CD1A|R7
913 CD1E - CD1A|R2
909 CD1A - CD1|FCB6|HTA1|R4|T6
912 CD1D - CD1A|R3
910 CD1B - CD1|CD1A|R1
9266 CYTH2 - ARNO|CTS18|CTS18.1|PSCD2|PSCD2L|SEC7L|Sec7p-L|Sec7p-like
30011 SH3KBP1 - CD2BP3|CIN85|GIG10|HSB-1|HSB1|MIG18
23607 CD2AP - CMS
89886 SLAMF9 - CD2F-10|CD2F10|CD84-H1|CD84H1|SF2001
10849 CD3EAP - ASE-1|ASE1|CAST|PAF49
445347 TARP - CD3G|TCRG|TCRGC1|TCRGC2
915 CD3D - CD3-DELTA|IMD19|T3D
920 CD4 - CD4mut
922 CD5L - AIM|API6|PRO229|SP-ALPHA|Spalpha
925 CD8A - CD8|Leu2|MAL|p32
927 CD8BP - CD8B2
54675 CRLS1 - C20orf155|CLS|CLS1|GCD10|dJ967N21.6
3681 ITGAD - ADB2|CD11D
3683 ITGAL - CD11A|LFA-1|LFA1A
3684 ITGAM - CD11B|CR3A|MAC-1|MAC1A|MO1A|SLEB6
3687 ITGAX - CD11C|SLEB6
290 ANPEP - APN|CD13|GP150|LAP1|P150|PEPN
115708 TRMT61A - C14orf172|GCD14|Gcd14p|TRM61|hTRM61
2526 FUT4 - CD15|ELFT|FCT3A|FUC-TIV|FUTIV|LeX|SSEA-1
2215 FCGR3B - CD16|CD16b|FCG3|FCGR3|FCR-10|FCRIII|FCRIIIb
4055 LTBR - CD18|D12S370|LT-BETA-R|TNF-R-III|TNFCR|TNFR-RP|TNFR2-RP|TNFR3|TNFRSF3
930 CD19 - B4|CVID3

 

 

 

 

 

 

 

30

自学CHIP-seq第二讲之过滤数据并比对

这个是有着非常成熟的流程了,我就不细讲了!

我们随机挑选两个文件来跑一下CHIP-seq的流程吧,其中一个是.部分进行免疫共沉淀前的DNA(input DNA)作为空白对照。

5.5G Apr 30 10:31 Xu_WT_rep2_BAF155.fastq

18G Feb 13 20:37 Xu_WT_rep2_Input.fastq

首先进行质量控制,过滤低质量的reads

这里我选取的是DynamicTrim.pl 和

脚本如下

for id in *fastq

do

echo $id

perl DynamicTrim.pl $id

done

接下来

for id in *.trimmed

do

echo $id

perl LengthSort.pl $id

Done

这样就得到了过滤后的reads,可以进行比对啦!

图片1

当然,中间文件可以删掉啦,不然太占空间了,我还只是取了两个数据,要是把这个文章的八个数据都跑完就太纠结了。

然后用bowtie比对

#samtools faidx hg19.fa

#Bowtie2-build hg19.fa hg19

for i in *single

do

bowtie2 -x /home/jmzeng/ref-database/hg19 -U $i -S  $i.sam

samtools view -bS $i.sam> $i.bam

done

输出的bam文件就需要用MASC这个软件来找peak了

30

自学CHIP-seq第一讲之文献解读

我这里选择的CHIP-seq文章题目是

CARM1 Methylates Chromatin Remodeling Factor BAF155 to Enhance Tumor Progression and Metastasis

文章链接http://www.sciencedirect.com/science/article/pii/S1535610813005369

这是2013年的文章,算是蛮新的了,主要探究了CARM1这个基因

然后我简单搜索了一些这个基因的信息

9606 10498 CARM1

- PRMT4

MIM:603934|HGNC:HGNC:23393|

Ensembl:ENSG00000142453|HPRD:09158|Vega:OTTHUMG00000180699

19 19p13.2 coactivator-associated arginine methyltransferase 1

protein-coding CARM1 coactivator-associated arginine methyltransferase histone-arginine methyltransferase CARM1|protein arginine N-methyltransferase 4 20150308

该基因是多种肿瘤相关的转录因子的共激活剂(激活蛋白;转录辅助激活蛋白;转录共同活化子)。

文章作者做了以下四件事

Knockout of CARM1 Using ZFN in Breast Cancer Cells

Identification of BAF155 as a Novel CARM1 Substrate

Methylation of BAF155 Promotes Tumor Growth and Metastasis

Methylated BAF155 Gains Unique Chromatin Association

 

所以就有两种细胞,一种是野生型WT,一种是突变的MUT细胞

然后它们分别做了两个重复,一种是input一种是BAF155免疫测序。

CHIP-seq一定是有一个input对照文件,和一个真正的免疫共沉淀的测序文件。

这样就有八个测序文件。

我们随机挑选两个文件来跑一下CHIP-seq的流程吧,其中一个是.部分进行免疫共沉淀前的DNA(input DNA)作为空白对照。

5.5G Apr 30 10:31 Xu_WT_rep2_BAF155.fastq

18G Feb 13 20:37 Xu_WT_rep2_Input.fastq

然后我随便在网上找了一个生信分析流程

  1. 标准信息分析
    a)   对测序数据进行base calling、raw data 数据整理及数据质量评估;
    b)   去接头污染,去低质量reads和产量情况统计
    c)   Bisulfite 测序序列与参考基因组序列的比对
    d)   深度和覆盖度分析
    e)   C 碱基的甲基化水平
    f)   全基因组甲基化水平分布趋势
  2. g)  全基因组DNA甲基化图谱
  3. h)  差异性甲基化区域(DMR)分析

 

参考

http://www.plob.org/2012/09/29/3760.html

http://www.plob.org/2012/01/09/1605.html

http://www.plob.org/2012/01/08/1538.html

 

30

阿里巴巴免费的服务器体验好差!

不知道为什么最近进入自己的网页后台总是很慢,发个日志也慢,很是郁闷!

本来以为是免费的空间快用完了,所以慢,结果一查,根本就没有用多,其实我很想投诉一下阿里巴巴!

想想该搞个国外服务器了,然后把网站搬家!

QQ截图20150430150101

30

Figtree的把进化树文件可视化

下载软件

http://tree.bio.ed.ac.uk/software/figtree/

我们这里就在window平台下使用,所以直接下载zip包即可,解压即可使用

准备数据

我这里就简单的用muscle生成了一个树文件来看看结果TRAV.fa 是一百多个类似的基因

muscle -in TRAV.fa -out tmp

muscle -maketree -in tmp  -out TRAV.tree

这个树文件TRAV.tree是Newick format,可以直接被figtree识别从而画图

软件使用

很简单,下载,点击即可使用,然后导入树文件,就可以直接出图啦!

Figtree的把进化树文件可视化368

30

Muscle进行多序列比对

软件的主页是

http://www.drive5.com/muscle/

进入主页,简单看看软件介绍,这个软件还是蛮牛的,一个人在家里自己写出来的,当然,对于普通人来说,这个软件跟clustalW没什么区别,反正都是多序列比对啦!

我们下载适合我们平台的版本即可!

Muscle进行多序列比对193

准备数据,我这里选择的是几个短小的蛋白

Muscle进行多序列比对215

 

这里有两种比对方式,都是很简单的命令

一种是先比对,再生成树文件(树的格式是Newick format, )

muscle -in mouse_J.pro -out mouse_J.pro.a

muscle -maketree -in mouse_J.pro.a -out mouse_J.phy (这里有两种构建树的方式)

另外一种是比对成aln格式的数据,然后用其它软件(phyml或者phylip)来生出树文件

muscle -in mouse_J.pro   -clwout seqs.aln

可以看到比对的效果还是蛮好的,是aln格式的比对文件,这个格式非常常用

Muscle进行多序列比对505

或者输出phy格式的比对文件,

muscle -in mouse_J.pro  -physout seqs.phy

Muscle进行多序列比对685

可以被phyml等软件识别,然后来构建进化树,见  http://www.bio-info-trainee.com/?p=626

21

美国Minnesota大学的生信全套课件分享

刚才在知乎什么看到了一篇分享pacbio的数据特征,顺便看到了Minnesota大学的关于生物信息的教程的ppt合集,所以就想打包下载。

https://www.msi.umn.edu/tutorial-materials

这个网页里面有64篇pdf格式的ppt,还有几个压缩包,本来是准备写爬虫来爬去的,但是后来想了想有点麻烦,而且还不一定会看,反正也是玩玩
就用linux的命令行简单实现了这个爬虫功能。
curl https://www.msi.umn.edu/tutorial-materials >tmp.txt
perl -alne '{/(https.*?pdf)/;print $1 if $1}' tmp.txt >pdf.address
perl -alne '{/(https.*?txt)/;print $1 if $1}' tmp.txt
perl -alne '{/(https.*?zip)/;print $1 if $1}' tmp.txt >zip.address
wget -i pdf.address
wget -i pdf.zip
这样就可以啦!
教程ppt列表如下,大家有兴趣的可以自行下载浏览。

2009-04-22-mrm-presentation_0.pdf               Matlab_viz_image_UMR.pdf
Analyzing ChIP at the command line.pdf          MaxQuant_Introduction_112409.pdf
Analyzing ChIP using Galaxy.pdf                 Maxquant-step-by-step_rs091124.pdf
Badalamenti_PacBio_tutorial_12-10-2014.pdf      MSI Applications Catalog Oct 21 MB slides.pdf
basics_chip_seq.pdf                             MSIIntro2013Jun18.pdf
Best_Practices_GATK_Variant_Detection_v1_0.pdf  MSIIntroBMEN5311.pdf
blast2go.pdf                                    MSI_Workshop_for_Introduction_to_Structure_based_Drug_Design.pdf
ClinProTools_0.pdf                              MTLB_GPUs.pdf
CUDA_Programming.pdf                            OpenMP.tutorial_1.pdf
cuda_tutorial_performance.pdf                   Open_Source_Proteomics_1.pdf
FLUENT_2009April21_final.pdf                    OptimizingWithGA.pdf
FLUENT_tutorial_2008aug14fin.pdf                Orbi_Data_Analysis_092811.pdf
galaxy_101_V4_ljm_0.pdf                         Partek Training Handout_miRNA and mRNA Data Analysis.pdf
GPU_tools.pdf                                   PerformanceTuning_itasca_11_27_12_0.pdf
gpututorial-msi.pdf                             PETSc_Tutorial.pdf
Hands_On_Tutorial_Using_ProTIP.pdf              Phi_Intro.pdf
Introduction to MSI Systems.pdf                 Protein_Grouping_FDR_Analysis_and_Database_Pratik_March2012_Draft.pdf
Introduction_to_PEAKS_0.pdf                     Proteomics_MSI_072309_Print.pdf
Introduction_to_SBDD.pdf                        pymol_v5.pdf
IntroMPI2011july19c.pdf                         QC_illumina_galaxy_V1_ljm.pdf
IntroMPI2012_July25-part1.pdf                   Quality Control of Illumina Data at the Command Line.pdf
IntroMSI2014.pdf                                remotevisualization.pdf
IntroNWChem.pdf                                 RISS_Hsapiens_variant_Detection_v3.0-small.pdf
IntroOpenMP_2011jun28b.pdf                      RNA_seq_Lecture2_2014_v2.pdf
Intro_to_GAMESS.pdf                             RNA-Seq mod1v6.pdf
IntroToGaussian09.pdf                           R_Spring2012_ver2.pdf
introtomolpro.pdf                               SchrodingerTutorial2011.pdf
Intro_to_MSI_Physicists.pdf                     Sybyl.pdf
intro-to-perl.pdf                               Tutorial-Hsap-v15.pdf
Matlab_11_29_UMR.pdf                            Tutorial-Stuber-v12-1.pdf
Matlab_PCT.pdf                                  unix2013.6.18.pdf
MATLAB_Tuning.pdf                               WRKSP_2_19.pdf

Total wall clock time: 40m 22s
Downloaded: 64 files, 249M in 40m 2s (106 KB/s)

我都已经下载好了,打包压缩到群里面啦!

21

HGNC数据库简介

人类基因命名委员会(HUGO Gene Nomenclature Committee);人类基因组命名委员会!

其实有了NCBI的entrez ID,然后还有refseq里面的ID,还有ensembl的ID,还有基因本身的功能英文缩略简称,已经很麻烦了,又来了一个HGNC,唉,头疼!

The HGNC approves both a short-form abbreviation known as a gene symbol, and also a longer and more descriptive name.

可以下载整个数据,用脚本慢慢研究研究

wget ftp://ftp.ebi.ac.uk/pub/databases/genenames/new/tsv/hgnc_complete_set.txt

 

还是看看BRCA1这个基因,里面的信息挺多的,主要看HGNC:1100,就是这个数据库对它这个基因的编号

HGNC:1100

BRCA1  这个是基因名,需要得到该组织的认可!!!!

breast cancer 1, early onset protein-coding gene gene with protein product Approved

17q21.31 17q21.31

"RNF53|BRCC1|PPP1R53|FANCS" "BRCA1/BRCA2-containing complex, subunit 1|protein phosphatase 1, regulatory subunit 53|Fanconi anemia, complementation group S" "Ring finger proteins|Protein phosphatase 1 regulatory subunits" "58|694"

1991-02-20T00:00:00Z

2015-04-18T00:00:00Z

672     这里是entrez ID

ENSG00000012048  这里是ensembl的ID,

OTTHUMG00000157426 uc002ict.3 U14680 NM_007294

"CCDS11453|CCDS11454|CCDS11455|CCDS11456|CCDS11459|CCDS11455|CCDS11456|CCDS11459|CCDS11454" P38398 1676470 MGI:104537 RGD:2218

"Breast Cancer|http://research.nhgri.nih.gov/bic/|BRCA1 database at LOVD-China|http://genomed.org/LOVD/BC/home.php?select_db=BRCA1|LOVD - Leiden Open Variation Database|http://chromium.liacs.nl/LOVD2/cancer/home.php?select_db=BRCA1|LOVD - Leiden Open Variation Database|http://proteomics.bio21.unimelb.edu.au/lovd/genes/BRCA1|LRG_292|http://www.lrg-sequence.org/LRG/LRG_292"

BRCA1 113705 119068

 

数据结构大概就是这个样子的了!

这几个数据库的内容都是互相链接的!

 

然后我们看看HGNC数据库的一些统计信息

http://www.genenames.org/cgi-bin/statistics

总共有40392个基因信息

其中18990个是能编码蛋白产物的基因,它们大多有GO注释

其中5927个是non-coding RNA,是现在的研究热门。

还有12546个是假基因,挺复杂的

最后还有1188个免疫相关基因,位置基因,病毒基因等等

 

最后,送给大家一个彩蛋!还有十一个物种也是有一个命名委员会的!

类似于 Mouse Gene Nomenclature Committee (MGNC).  Please see the following links:

 

参考文献;

Gray KA, Yates B, Seal RL, Wright MW, Bruford EA. genenames.org: the HGNC resources in 2015. Nucleic Acids Res. 2015 Jan;43(Database issue):D1079-85. doi: 10.1093/nar/gku1071. PMID:25361968

18

Bowtie算法第六讲-tally法对bwt索引进行搜索

因为要讲搜索,所以我选择了一个长一点的字符串来演示多种情况的搜索

perl rotation_one_by_one.pl atgtgtcgtagctcgtnncgt

程序运行的结果如下

$ATGTGTCGTAGCTCGTNNCGT 21

AGCTCGTNNCGT$ATGTGTCGT 9

ATGTGTCGTAGCTCGTNNCGT$ 0

CGT$ATGTGTCGTAGCTCGTNN 18

CGTAGCTCGTNNCGT$ATGTGT 6

CGTNNCGT$ATGTGTCGTAGCT 13

CTCGTNNCGT$ATGTGTCGTAG 11

GCTCGTNNCGT$ATGTGTCGTA 10

GT$ATGTGTCGTAGCTCGTNNC 19

GTAGCTCGTNNCGT$ATGTGTC 7

GTCGTAGCTCGTNNCGT$ATGT 4

GTGTCGTAGCTCGTNNCGT$AT 2

GTNNCGT$ATGTGTCGTAGCTC 14

NCGT$ATGTGTCGTAGCTCGTN 17

NNCGT$ATGTGTCGTAGCTCGT 16

T$ATGTGTCGTAGCTCGTNNCG 20

TAGCTCGTNNCGT$ATGTGTCG 8

TCGTAGCTCGTNNCGT$ATGTG 5

TCGTNNCGT$ATGTGTCGTAGC 12

TGTCGTAGCTCGTNNCGT$ATG 3

TGTGTCGTAGCTCGTNNCGT$A 1

TNNCGT$ATGTGTCGTAGCTCG 15

它的BWT及索引是

T 21

T 9

$ 0

N 18

T 6

T 13

G 11

A 10

C 19

C 7

T 4

T 2

C 14

N 17

T 16

G 20

G 8

G 5

C 12

G 3

A 1

G 15

然后得到它的tally文件如下

图片1

接下来用我们的perl程序在里面找字符串

第一次我测试 GTGTCG 这个字符串,程序可以很清楚的看到它的查找过程。

perl search_char.pl    GTGTCG   tm.tally

your last char is G

start is 7 ; and end is 13

now it is number 5 and the char is C

start is 3 ; and end is 6

now it is number 4 and the char is T

start is 17 ; and end is 19

now it is number 3 and the char is G

start is 10 ; and end is 11

now it is number 2 and the char is T

start is 19 ; and end is 20

now it is number 1 and the char is G

start is 11 ; and end is 12

It is just one perfect match !

The index is 2

第二次我测试一个多重匹配的字符串GT,在原字符串出现了五次的

perl search_char.pl  GT  tm.tally

your last char is T

start is 15 ; and end is 22

now it is number 1 and the char is G

start is 8 ; and end is 13

we find more than one perfect match!!!

8 13

One of the index is 11

One of the index is 10

One of the index is 19

One of the index is 7

One of the index is 4

One of the index is 2

One of the index is 14

惨了,这个是很严重的bug,不知道为什么,对于多个匹配总是会多出那么一点点的结果。

去转换矩阵里面查看,可知,前面两个结果11和10是错误的。

CTCGTNNCGT$ATGTGTCGTAG 11

GCTCGTNNCGT$ATGTGTCGTA 10

GT$ATGTGTCGTAGCTCGTNNC 19

GTAGCTCGTNNCGT$ATGTGTC 7

GTCGTAGCTCGTNNCGT$ATGT 4

GTGTCGTAGCTCGTNNCGT$AT 2

GTNNCGT$ATGTGTCGTAGCTC 14

最后我们测试未知字符串的查找。

perl search_char.pl ACATGTGT tm.tally

your last char is T

start is 15 ; and end is 22

now it is number 7 and the char is G

start is 8 ; and end is 13

now it is number 6 and the char is T

start is 19 ; and end is 21

now it is number 5 and the char is G

start is 11 ; and end is 12

now it is number 4 and the char is T

start is 20 ; and end is 21

now it is number 3 and the char is A

start is 2 ; and end is 3

now it is number 2 and the char is C

start is 3 ; and end is 3

we can just find the last 6 char ,and it is ATGTGT

原始字符串是ATGTGTCGTAGCTCGTNNCGT,所以查找的挺对的!!!

 

[perl]

$a=$ARGV[0];

$a=uc $a;

open FH,"<$ARGV[1]";

while(<FH>){

chomp;

@F=split;

$hash_count_atcg{$F[0]}++;

$hash{$.}=$_;

# the first line is $ and the last char and the last index !

}

$all_a=$hash_count_atcg{'A'};

$all_c=$hash_count_atcg{'C'};

$all_g=$hash_count_atcg{'G'};

$all_n=$hash_count_atcg{'N'};

$all_t=$hash_count_atcg{'T'};

#print "$all_a\t$all_c\t$all_g\t$all_t\n";

$len_a=length $a;

$end_a=$len_a-1;

#print "your query is $a\n";

#print "and the length of your query is $len_a \n";

$after=substr($a,$end_a,1);

#we fill search your query from the last char !

if ($after eq 'A') {

$start=2;

$end=$all_a+1;

}

elsif ($after eq 'C') {

$start=$all_a+1;

$end=$all_a+$all_c+1;

}

elsif ($after eq 'G') {

$start=$all_a+$all_c+1;

$end=$all_a+$all_c+$all_g+1;

}

elsif ($after eq 'T'){

$start=$all_a+$all_c+$all_g+$all_n+1;

$end=$all_a+$all_c+$all_g+$all_t+$all_n+1;

}

else {die "error !!! we just need A T C G !!!\n"}

print "your last char is $after\n ";

print "start is $start ; and end is $end \n";

foreach (reverse (1..$end_a)){

$after=substr($a,$_,1);

$before=substr($a,$_-1,1);

($start,$end)=&find_level($after,$before,$start,$end);

print "now it is number $_ and the char is $before \n ";

print "start is $start ; and end is $end \n";

if ($_  > 1 && $start == $end) {

$find_char=substr($a,$_);

$find_len=length $find_char;

print "we can just find the last $find_len char ,and it is $find_char \n";

#return "miss";

last;

}

if ($_ == 1) {

if (($end-$start)==1) {

print "It is just one perfect match ! \n";

my @F_start=split/\s+/,$hash{$end};

print "The index is $F_start[1]\n";

#return $F_start[1];

last;

}

else {

print "we find more than one perfect match!!!\n";

print "$start\t$end\n";

foreach  (($start-1)..$end) {

my @F_start=split/\s+/,$hash{$_};

print "One of the index is $F_start[1]\n";

}

#return "multiple";

last;

}

}

}

sub find_level{

my($after,$before,$start,$end)=@_;

my @F_start=split/\s+/,$hash{$start};

my @F_end=split/\s+/,$hash{$end};

if ($before eq 'A') {

return ($F_start[2]+1,$F_end[2]+1);

}

elsif ($before eq 'C') {

return ($all_a+$F_start[3]+1,$all_a+$F_end[3]+1);

}

elsif ($before eq 'G') {

return ($all_a+$all_c+1+$F_start[4],$all_a+$all_c+1+$F_end[4]);

}

elsif ($before eq 'T') {

return ($all_a+$all_c+$all_g+$all_n+1+$F_start[5],$all_a+$all_c+$all_g+1+$all_n+$F_end[5]);

}

else {die "error !!! we just need A T C G !!!\n"}

}

[/perl]

 

原始字符串是atgtgtcgtagctcgtnncgt

 

18

Bowtie算法第五讲-index2tally

前面讲到了如何用笨方法进行字符串搜索,也讲了如何构建bwt索引,和把bwt索引还原成字符串!

原始字符串是ATGCGTANNGTC

排序过程是下面的

$ATGCGTANNGTC 12

ANNGTC$ATGCGT 6

ATGCGTANNGTC$ 0

C$ATGCGTANNGT 11

CGTANNGTC$ATG 3

GCGTANNGTC$AT 2

GTANNGTC$ATGC 4

GTC$ATGCGTANN 9

NGTC$ATGCGTAN 8

NNGTC$ATGCGTA 7

TANNGTC$ATGCG 5

TC$ATGCGTANNG 10

TGCGTANNGTC$A 1

现在讲讲如何根据bwt索引构建tally,并且用tally搜索方法来搜索字符串!

首先是bwt索引转换为tally

C 12

T 6

$ 0

T 11

G 3

T 2

C 4

N 9

N 8

A 7

G 5

G 10

A 1

这个其实非常简单的,tally就是增加四列计数的列即可

[perl]

$hash_count{'A'}=0;

$hash_count{'C'}=0;

$hash_count{'G'}=0;

$hash_count{'T'}=0;

open FH ,"<$ARGV[0]";

while(<FH>){

        chomp;

@F=split;

$last=$F[0]; # 读取上面的tally文件,分列,判断第一列,并计数

        $hash_count{$last}++;

   print  "$_\t$hash_count{'A'}\t$hash_count{'C'}\t$hash_count{'G'}\t$hash_count{'T'}\n";

}

[/perl]

输出的tally如下

C 12 0 1 0 0

T 6 0 1 0 1

$ 0 0 1 0 1

T 11 0 1 0 2

G 3 0 1 1 2

T 2 0 1 1 3

C 4 0 2 1 3

N 9 0 2 1 3

N 8 0 2 1 3

A 7 1 2 1 3

G 5 1 2 2 3

G 10 1 2 3 3

A 1 2 2 3 3

接下来就是针对这个tally的查询函数了

 

18

Bowtie 算法第四讲

由于之前就简单的看了看bowtie作者的ppt,没有完全吃透就开始敲代码了,写了十几个程序最后我自己都搞不清楚进展到哪一步了,所以我现在整理一下,从新开始!!!

 

首先,bowtie的作用就是在一个大字符串里面搜索一个小字符串!那么本身就有一个非常笨的复杂方法来搜索,比如,大字符串长度为100万,小字符串为10,那么就依次取出大字符串的10个字符来跟小字符串比较即可,这样的算法是非常不经济的,我简单用perl代码实现一下。

[perl]

#首先读取大字符串的fasta文件

open FH ,"<$ARGV[0]";

$i=0;

while (<FH>) {

next if /^>/;

chomp;

$a.=(uc);

}

#print "$a\n";

#然后接受我们的小的查询字符串

$query=uc $ARGV[1];

$len=length $a;

$len_query=length $query;

$a=$a.'$'.$a;

#然后依次循环取大字符串来精确比较!

foreach (0..$len-1){

if (substr($a,$_,$len_query) eq $query){

print "$_\n";

#last;

}

}

[/perl]

 

这样在时间复杂度非常恐怖,尤其是对人的30亿碱基。

 

正是因为这样的查询效率非常低,所以我们才需要用bwt算法来构建索引,然后根据tally来进行查询

其中构建索引有三种方式,我首先讲最效率最低的那种索引构造算法,就是依次取字符串进行旋转,然后排序即可。

[perl]

$a=uc $ARGV[0];

$len=length $a;

$a=$a.'$'.$a;

foreach (0..$len){

$hash{substr($a,$_,$len+1)}=$_;

}

#print "$_\t$hash{$_}\n" foreach sort keys %hash;

print  substr($_,-1),"\t$hash{$_}\n" foreach sort keys %hash;

[/perl]

这个算法从时间复杂度来讲是非常经济的,对小字符串都是瞬间搞定!!!

perl rotation_one_by_one.pl atgcgtanngtc 这个字符串的BWT矩阵索引如下!

C 12

T 6

$ 0

T 11

G 3

T 2

C 4

N 9

N 8

A 7

G 5

G 10

A 1

但同样的,它也有一个无法避免的弊端,就是内存消耗太恐怖。对于30亿的人类碱基来说,这样旋转会生成30亿乘以30亿的大矩阵,一般的服务器根本hold不住的。

 

最后我讲一下,这个BWT矩阵索引如何还原成原字符串,这个没有算法的差别,因为就是很简单的原理。

[perl]

#first read the tally !!!

#首先读取上面输出的BWT矩阵索引文件。

open FH,"<$ARGV[0]";

$hash_count{'A'}=0;

$hash_count{'C'}=0;

$hash_count{'G'}=0;

$hash_count{'T'}=0;

while(<FH>){

        chomp;

        @F=split;

        $hash_count{$F[0]}++;

        $hash{$.}="$F[0]\t$F[1]\t$hash_count{$F[0]}";

#print "$hash{$.}\n";

}

$all_a=$hash_count{'A'};        

$all_c=$hash_count{'C'};        

$all_g=$hash_count{'G'};        

$all_t=$hash_count{'T'};

$all_n=$hash_count{'N'};

#start from the first char !

$raw='';

&restore(1);

sub restore{

my($num)=@_;

my @F=split/\t/,$hash{$num};

$raw.=$F[0];

   my $before=$F[0];

     if ($before eq 'A') { 

$new=$F[2]+1;

        }

        elsif ($before eq 'C') {

               $new=1+$all_a+$F[2];

        }

        elsif ($before eq 'G') {

               $new=1+$all_a+$all_c+$F[2];

        }

elsif ($before eq 'N') {

                $new =1+$all_a+$all_c+$all_g+$F[2];

        }

        elsif ($before eq 'T') {

                $new=1+$all_a+$all_c+$all_g+$all_n+$F[2];

        }

        elsif ($before eq '$') {

chop $raw;

                $raw = reverse $raw;

print "$raw\n";

exit;

        }

else {die "error !!! we just need A T C N G !!!\n"}

#print "$F[0]\t$new\n";

&restore($new);

}

[/perl]

 

 

16

Annovar使用记录

至于如何安装该软件,请见上一个教程

一.首先把snp-calling步骤的VCF文件转为annovar软件要求的格式

convert2annovar.pl   -format vcf4   12.vcf >12.annovar

Annovar使用记录108

二.进行注释

命令行参数比较多,还是用脚本来运行

# define path

infolder=/home/jmzeng/hoston/diff

outfolder=$infolder

annovardb=/home/jmzeng/bio-soft/annovar/humandb

# start annotating

/home/jmzeng/bio-soft/annovar/annotate_variation.pl \

--buildver hg19 \

--geneanno \

--outfile ${outfolder}/12.anno \

${infolder}/12.annovar  \

${annovardb}

三.输出结果解读

2.6M Apr 14 22:32 12.anno.exonic_variant_function

1.9K Apr 14 22:32 12.anno.log

1.3M Apr 14 22:32 12.anno.variant_function

重点是后缀为exonic_variant_function,这个文件对每一个vcf的突变都进行了注释。

Annovar使用记录617

这个结果就可以用来解析了,可以根据实验设计来找到自己感兴趣的突变。

第5.6列是染色体及pos坐标

第4列信息非常复杂,是突变的注释

第12列是测序深度,一般要大于20

我这里是先把注释文件转换成以下格式

location:chr1:874467 SAMD11:NM_152486:exon6:c.G478A:p.D160N

location:chr1:888639 NOC2L:NM_015658:exon9:c.A918G:p.E306E

location:chr1:888659 NOC2L:NM_015658:exon9:c.A898G:p.I300V

location:chr1:916549 PERM1:NM_001291367:exon2:c.T58C:p.W20R

location:chr1:949608 ISG15:NM_005101:exon2:c.G248A:p.S83N

location:chr1:980552 AGRN:NM_198576:exon13:c.G2266A:p.A756T

location:chr1:1114699 TTLL10:NM_001130045:exon4:c.G104A:p.R35Q

location:chr1:1158631 SDF4:NM_016176:exon4:c.T570C:p.D190D

location:chr1:1158631 SDF4:NM_016547:exon4:c.T570C:p.D190D

location:chr1:1164073 SDF4:NM_016176:exon2:c.C101T:p.A34V

然后比较两个文件,取不同的突变来格式化输出。

16

推荐linux学习博客-每日一linux命令

竹子-博客(.NET/Java/Linux/架构/管理/敏捷)

思索、感悟、践行!走向高效,快乐,平衡!

http://www.cnblogs.com/peida/tag/%E6%AF%8F%E6%97%A5%E4%B8%80linux%E5%91%BD%E4%BB%A4/default.html?page=1

已下目录是本人用爬虫爬取的!

每天一个linux命令(61):wget命令

每天一个linux命令(60):scp命令

每天一个linux命令(59):rcp命令

每天一个linux命令(58):telnet命令

每天一个linux命令(57):ss命令

每天一个linux命令(56):netstat命令

每天一个linux命令(55):traceroute命令

每天一个linux命令(54):ping命令

每天一个linux命令(53):route命令

每天一个linux命令(52):ifconfig命令

每天一个linux命令(51):lsof命令

每天一个linux命令(50):crontab命令

每天一个linux命令(49):at命令

每天一个linux命令(48):watch命令

每天一个linux命令(47):iostat命令

每天一个linux命令(46):vmstat命令

每天一个linux命令(45):free 命令

每天一个linux命令(44):top命令

每天一个linux命令(43):killall命令

每天一个linux命令(42):kill命令

每天一个linux命令(41):ps命令

每天一个linux命令(40):wc命令

每天一个linux命令(39):grep 命令

每天一个linux命令(38):cal 命令

每天一个linux命令(37):date命令

每天一个linux命令(36):diff 命令

每天一个linux命令(35):ln 命令

每天一个linux命令(34):du 命令

每天一个linux命令(33):df 命令

每天一个linux命令(32):gzip命令

每天一个linux命令(31): /etc/group文件详解

每天一个linux命令(30): chown命令

每天一个linux命令(29):chgrp命令

每天一个linux命令(28):tar命令

每天一个linux命令(27):linux chmod命令

每天一个linux命令(26):用SecureCRT来上传和下载文件

每天一个linux命令(25):linux文件属性详解

每天一个linux命令(24):Linux文件类型与扩展名

每天一个linux命令(23):Linux 目录结构

每天一个linux命令(22):find 命令的参数详解

每天一个linux命令(21):find命令之xargs

每天一个linux命令(20):find命令之exec

每天一个linux命令(19):find 命令概览

每天一个linux命令(18):locate 命令

每天一个linux命令(17):whereis 命令

每天一个linux命令(16):which命令

每天一个linux命令(15):tail 命令

每天一个linux命令(14):head 命令

每天一个linux命令(13):less 命令

每天一个linux命令(12):more命令

每天一个linux命令(11):nl命令

每天一个linux命令(10):cat 命令

每天一个linux命令(9):touch 命令

每天一个linux命令(8):cp 命令

每天一个linux命令(7):mv命令

每天一个linux命令(6):rmdir 命令

每天一个linux命令(5):rm 命令

每天一个linux命令(4):mkdir命令

每天一个linux命令(3):pwd命令

每天一个linux命令(2):cd命令

每天一个linux命令(1):ls命令

 

15

human已经被研究的snp竟然有一亿多个?

我在NCBI里面下载了一个dbsnp_142数据库文件,发现它居然有2.5G的大小,我感到很不可思议,毕竟人的基因组也就3G,就30亿的碱基嘛。研究过的突然竟然有110,917,213 ,高达一亿个!!!

谁能给我解释一下呢!

而且人只有十万多个蛋白,2.2万多个基因!

jmzeng@ubuntu:/home/jmzeng/hoston/diff/snp$ wc -l dbsnp_142_chrom_id_rs
110917213 dbsnp_142_chrom_id_rs
jmzeng@ubuntu:/home/jmzeng/hoston/diff/snp$ tail dbsnp_142_chrom_id_rs
MT    16429    rs150751410
MT    16443    rs371960162
MT    16456    rs142662828
MT    16482    rs386419986
MT    16497    rs376846509
MT    16512    rs373943637
MT    16519    rs3937033
MT    16526    rs386829315
MT    16527    rs386829316
MT    16529    rs370705831
jmzeng@ubuntu:/home/jmzeng/hoston/diff/snp$ head dbsnp_142_chrom_id_rs
1    10108    rs62651026
1    10109    rs376007522
1    10139    rs368469931
1    10144    rs144773400
1    10150    rs371194064
1    10177    rs201752861
1    10177    rs367896724
1    10180    rs201694901
1    10228    rs143255646
1    10228    rs200462216

15

Vcf文件的突变ID号注释

VCF是1000genome计划定义的测序比对突变说明文件,熟悉VCF文件的都知道,第三列是ID号,也就是说对该突变在dbsnp的数据库的编号。大多时候都是用点号占位,代表不知道在dbsnp的数据库的编号,这时候就需要我们自己来注释了。

Vcf文件的突变ID号注释134

其实,这是一个非常简单的事情,因为有了CHROM和pos,只要找到一个文件,就可以自己写脚本来映射到它的ID号,但是找这个文件比较困难,我也是搜索了好久才找到的。

http://varianttools.sourceforge.net/Annotation/DbSNP

这里面提到了最新版的数据库是dbSNP138

The default version of our dbSNP annotation is currently referring to dbSNP138 (using hg19 coordinates) as shown below. However, users can also retrieve older versions of dbSNP: db135, dbSNP129, dbSNP130, dbSNP131 and dbSNP132. The 129 and 130 versions use hg18 as a reference genome and 131, 132, 135 and later use hg19. The archived versions can be used by a variant tools project by referring to their specific names - for example: dbSNP-hg18_129.

所以我就换了关键词,终于搜的了

http://www.ncbi.nlm.nih.gov/projects/SNP/snp_summary.cgi?view+summary=view+summary&build_id=138

http://asia.ensembl.org/info/genome/variation/sources_documentation.html?redirect=no

SNP 138 database (232,952,851 million altogether).

Vcf文件的突变ID号注释1276

有一个bioconductor包是专门来做snp过滤的

http://www.bioconductor.org/packages/release/bioc/html/VariantAnnotation.html

首先下载vcf文件。

nohup wget ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/00-All.vcf.gz &

这个文件很大,解压开是

如果大家对snp不了解,可以去查看它的各种介绍以及分类

http://moma.ki.au.dk/genome-mirror/cgi-bin/hgTrackUi?db=hg19&g=snp138

 

其实我这里本来有个hg19_snp141.txt文件,如下

1 10020 A - .

1 10108 C T .

1 10109 A T .

1 10139 A T .

1 10145 A - .

1 10147 C - .

1 10150 C T .

1 10177 A C .

1 10180 T C .

1 10229 A - .

 

还可以下载一些文件,如bed_chr_1.bed

chr1 175292542 175292543 rs171 0 -

chr1 20542967 20542968 rs242 0 +

chr1 6100897 6100898 rs538 0 -

chr1 93151988 93151989 rs546 0 +

chr1 15220328 15220329 rs549 0 +

chr1 203744004 203744005 rs568 0 +

chr1 23854550 23854551 rs665 0 -

chr1 53213656 53213657 rs672 0 +

chr1 173907422 173907423 rs677 0 -

当然还有那个ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/00-All.vcf.gz  18G的文件,是VCF格式

##fileformat=VCFv4.0

##fileDate=20150218

##source=dbSNP

##dbSNP_BUILD_ID=142

##reference=GRCh38

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO

1       10108   rs62651026      C       T       .       .       RS=62651026;RSPOS=10108;dbSNPBuildID=129;SSR=0;SAO=0;VP=0x050000020005000002000100;WGT=1;VC=SNV;R5;ASP

所以这个文件就是我们想要的最佳文件,提取前三列就够啦

#CHROM  POS     ID

1 10108 rs62651026

1 10109 rs376007522

1 10139 rs368469931

1 10144 rs144773400

1 10150 rs371194064

1 10177 rs201752861

1 10177 rs367896724

1 10180 rs201694901

1 10228 rs143255646

1 10228 rs200462216

这样就可以通过脚本用hash把我们自己找到的hash跟数据库的rs编号对应起来啦

 

13

用phyML对多重比对phy文件来构建进化树

本来还没这么快写进化专题的,但是有个朋友要我帮忙跑一下她的phy文件,看看能否生成树,我就试用了一下phyml这个软件,挺简单的。

一、下载并安装该软件

这是一个很简单的软件,我们直接下载它的二进制程序就可以直接使用啦,官网里面的压缩包里面有各种平台的二进制程序,我这里用linux64的。

wget http://www.atgc-montpellier.fr/download/binaries/phyml/PhyML-3.1.zip

unzip PhyML-3.1.zip

构建进化树phyml241

二.准备文件

它需要的phy格式的多重比对结果文件,一般是clustalW或者muscle比对的结果

构建进化树phyml294

可以看到是53个蛋白,多重比对后的公共序列长度是325个氨基酸。

三.命令

./PhyML-3.1_linux64   -i test.phy  -d aa   -b 1000   -m LG   -f m -v e -a e -o tlr

这些参数在运行的时候都会显示出来

构建进化树phyml434

具体解释见博客 http://www.chenlianfu.com/?p=2221

 

四,输出文件

这个时间会很久,大家有心里准备!!!总共会输出四个文件,

test.phy_phyml_tree.txt        :    最大似然法构建的进化树

test.phy_phyml_boot_stats.txt  :    bootstrap 的统计信息

test.phy_phyml_boot_trees.txt  :    bootstrap 树

test.phy_phyml_stats.txt       :    程序运行的中的参数和结果统计

 

然后我们的那个test.phy_phyml_tree.txt  就可以用figtree等软件画图啦!!!

 

 

 

 

 

 

 

 

 

10

查找某个基因上面的snp位点

进入网页 http://www.ncbi.nlm.nih.gov/projects/SNP/

image001

其实就是http://www.ncbi.nlm.nih.gov/snp 这个网页

image003

可以看到这个基因上面发表的snp非常多,共有14893个。

每个突变的各种信息都很完全,比如第一个snp位点, 它的ID是rs12516,在BRCA1基因上面。是17号染色体的43044391的碱基突变,在refseq数据库里面显示有两个NG,一个NC,还有五个NM都突变了,还有一堆XM就无所谓啦。

http://www.ncbi.nlm.nih.gov/projects/SNP/snp_ref.cgi?rs=12516

image005

而且是有文献支持的,在1000genomes计划里面也有发表。而且是hg19和hg38里面是不同的坐标!

发表它的文献是 Associations between single nucleotide polymorphisms in double-stranded DNA repair pathway genes and familial breast cancer.

 

 

10

一步一步运行软件系列合集

这些是很久以前写的一些教程,是关于进化树构建和全基因组关联分析的!

gwas-plink分析教程.pdf
plink的统计基础.ppt
一步一步构建系统进化树.pdf
一步一步运行blast.pdf
一步一步运行inparanoid蛋白聚类.pdf
一步一步运行PLINK-part1.pdf
一步一步运行plink-part2.pdf
用PhyML构建系统发育树.pptx
进化树的构建分子原理.pdf

都在云盘(http://pan.baidu.com/s/1jIvwRD8 )里面,群空间(201161227)里面也有!

暂时应该不会写这些教程了,因为没有项目,实在没有动力去做那么多事情