CENTOS 6.0 开放端口80

CENTOS 6.0 开放端口80 3306 22等 iptables

博客分类： centos
iptablescentos
#/sbin/iptables -I INPUT -p tcp --dport 80 -j ACCEPT
#/sbin/iptables -I INPUT -p tcp --dport 22 -j ACCEPT
#/sbin/iptables -I INPUT -p tcp --dport 3306 -j ACCEPT
然后保存：
#/etc/rc.d/init.d/iptables save

查看打开的端口：
# /etc/init.d/iptables status

-------------------------------------------------------
补充说明：

#关闭防火墙
/etc/init.d/iptables stop
service iptables stop # 停止服务
#查看防火墙信息
/etc/init.d/iptables status

#开放端口:8080
/sbin/iptables -I INPUT -p tcp --dport 8080 -j ACCEPT
#重启防火墙以便改动生效:(或者直接重启系统)
/etc/init.d/iptables restart
#将更改进行保存
/etc/rc.d/init.d/iptables save

另外直接在/etc/sysconfig/iptables中增加一行：
-A RH-Firewall-1-INPUT -m state –state NEW -m tcp -p tcp –dport 8080 -j ACCEPT
#永久关闭防火墙
chkconfig –level 35 iptables off #此方法源自网络，未实验，安全考虑拒绝使用此方法

Rocks clusters安装经验总结

Rocks clusters安装经验总结

2009-08-12 23:03:09| 分类： 技术 |字号 订阅
经过几天的苦苦摸索，终于搞定了Rocks clusters的安装和初步配置。现将安装过程中遇到的问题及解决经验做一总结，以便给自己以后重新安装时做一参考，或者能帮助也想安装Rocks cluster的网友顺利一些，能节省一些花在摸索问题解决办法上的解决时间。

很早一段时间以前就下载了oscar集群安装软件，虽然很仔细的阅读了安装指南，但在安装的过程中还是遇到了很多问题，而不能顺利的成功安装，也花时间在网上找了很多零碎的文档，最后还是没能安装成功。在网上找文档时，得知了Rocks clusters也是一个非常流行的集群软件，还是免费开放的，并能实现自动安装，免去了修改很多配置文件的烦恼（在以前搭建PVM集群时就需要修改挺多配置文件，想起就有点头大），于是到其官方网站http://www.rocksclusters.org下载了当时最新版本5.1的Rocks安装程序（现在最新版本是5.2），我用来搭建集群的机器都是32位的i386类型的，于是下载了适用于i386的5张安装光盘，1张DVD光盘和4张CD光盘，DVD光盘上的东西（包括kernerl、base、os-1、os-2）都是必装的，其他几张CD上的(os-3到os-6)是选装的。在安装指南提到，Frontend节点机必须要有2块网卡才行，其中一个网卡接其他集群节点机，另一个网卡接外网。但我安装集群的目的是自己来做并行计算用，不打算与外网连通，所以就想当然的认为Frontend节点机有一个网卡就可以了，但正是因为这个想当然，使我在后面的安装过程中遇到了很多问题，浪费了很多时间去摸索，当意识到必须有两块网卡才可以之后，原来的很多问题也都自动消失了。所以，Frontend节点机必须有两块网卡才能安装Rocks成功，切记切记！

其实在我意识到Frontend节点机必须要有两块网卡之前，我也已经顺利安装了Frontend节点机和一台compute节点机，当时这样安装的Rocks系统是存在挺多问题的，比如qmaster等服务无法启动、也无法采用VM虚拟机的方式来作为计算节点。这些也是很长一点时间以后在重新来接着安装Rocks集群时才知道的。

安装Rocks cluster的步骤如下：

1、首先需要安装Frontend节点机，相当于一个集群的服务节点机：

（1）在CMOS中将光驱设为开机启动第一设备，将包含kernel包的Rocks集群的 DVD安装盘或只包含kernel包的CD光盘放入光驱，开机启动。

（2）当出现一个Rocks安装boot:界面时，会提示输入build加回车来开始安装Frontend节点（如果在boot:提示符后直接回车的话是开始安装compute或vm container等节点机，在后面安装compute或vm container节点机时会用到）。

（3）设置Frontend节点机中连接外网的网卡eth1的IP地址等网络配置。如果连接外网的网卡是采用DHCP的方式来获取IP的话，就可以跳过这步，系统会给eth1网卡自动分配IP并设置网络。如果采用固定IP的话，就需要手动设置IP地址及掩码、DNS等。

（4）接下来会提示选择安装文件的来源（CD/DVD或官方网址），我选择的是CD/DVD方式，然后开始选择需要安装的文件包（rocks中称为roll)，这时需要将DVD光盘重新放入或kenerl、base等、os-1，is-2这4张CD光盘依次放入，把需要安装的roll都打勾选择安装，这些roll都是必须要选的，如果没有选择全，在安装的过程中就会出现错误，只能重新启动安装。当然这时候也可以选择os-3至os-4这几张光盘里的roll来安装。等选择完需要安装的roll以后，就可以点next来继续安装了。

（5）这时需要选择对硬盘的分区方式（自动或手动），选择“自动”就会以默认的分区格式来对整块硬盘进行分区，选择“手动”就会进入分区工具界面，具体的分区要求可参照rocks安装说明上的推荐。我选择的是“手动”，因为“自动”方式会默认rocks安装占用整块硬盘，而我的硬盘上还装有WindowsXP系统，所以采用“手动”方式使得Rocks安装只占用其中一部分硬盘，这样Frontend节点机中也可以有双系统了。分区格式化硬盘以后就开始进入下面的安装过程了。

（6）提示依次放入包含前面选择安装的roll的各张光盘，会自动将各种光盘上的内容都拷贝到硬盘上，然后出现安装进度条开始安装过程。等一段时间安装完成后，Frontend节点机的安装就算大功告成了！！！

2. 然后是其他节点机（compute或vm container等）的安装

（1）以root登录Frontend节点机，在命令行运行insert-ethers，这时会弹出一个界面让选择是安装compute节点还是vm container节点或其他类型的节点。选择了要安装的节点机类型后，点击OK后，就会进入等待节点机安装请求的界面。这时就可以启动要安装的节点机开始安装过程了。

（2）安装时可以用光盘启动的方式或PXE网络启动的方式启动节点机，PXE网络启动方式需要在CMOS里设置，由于我的机器主板不支持PXE启动方式，所以我选择的是光盘启动，将包含kernel包的DVD或CD光盘放入要安装的节点机的光驱，引导启动节点机。

（3）这时会进入和安装Fontend节点机时同样的boot:提示符界面，现在不需要敲入build，而是直接回车就可以了，接下来就会自动进行整个安装过程，包括向Frontend节点机请求IP地址、从Frontend节点机拷贝kickstart安装文件及各种roll文件等。可以在Frontend节点机上运行"rocks-console compute-0-0(正在安装的节点机名) “来监控安装过程。

（4）节点机安装完就会自动重启，重启完成后可以在Frontend节点机运行rocks list host命令来列出Rocks cluster中已经安装好的节点机，如果能看到刚才安装的节点机就说明已经正确安装。

（5）重复以上的步骤可以安装其他的节点机。

在安装过程中我发现了Rocks cluster提供了通过VM虚拟机来构建集群的方法，这样我就可以只用两台机器来模拟包含多个节点机的集群了，从而方便在这个虚拟集群上来开发调试并行程序，等并行程序开发好以后再转移到包含多个物理节点机的集群上来测试并行计算的性能和速度。

Rocks cluster通xen技术来实现虚拟机，首先需要采用前面介绍的步骤来安装Frontend节点机，要保证在安装过程中选择安装了xen包，其他过程和前面都相同。然后再以和安装compute节点机类似的步骤安装vm container节点机，安装完以后可以用rocks list host 命令来查看是否安装成功，一般vm container节点机的主机名称会自动命名为"vm-container-0-0"这种形式。

虚拟集群的方法根据Frontend节点机是否虚拟机可分为两种，一种是在Frontend节点机也是作为虚拟机来安装，一种是Frontend节点机和前面的相同，都是一般的物理节点机。这两种方法在xen roll安装指南中都有详细介绍，我选用的是后一种方法。需要提一下的是，我安装的是最新的版本5.2，但在一开始安装虚拟集群时遇到了一些问题，不能成功安装compute虚拟机节点，后来在网上Rocks讨论区里查了一些文档才发现5.2版本在这方面存在bug，并且开发者已经写好了修正的fix包，名称为service-pack 5.2.1，官方网站上已经可以下载，在安装Frontend节点机时，这个包也要一起安装上才行。

我采用的后一种方式的Rocks虚拟集群的具体安装步骤是：

（1）在安装完Frontend节点机和vm container节点机后，在Frontend节点机上运行：

rocks add host vm vm-container-0-0 membership="Compute"

来在节点机vm-contianer-0-0上安装一个VM虚拟机作为compute类型的节点机，新安装的虚拟机的主机名一般会自动命名为“compute-0-0-0"这种形式，当然这个命令还可以加入其它一些选项，比如虚拟机的名称、分配内存大小等等，具体可参见xen roll安装指南。

（2）然后在Frontend节点机上运行：

rocks start host vm compute-0-0-0 install=yes

来启动刚才安装的虚拟机compute-0-0-0，并开始compute节点机的安装过程。install=yes一定要加上，否则无法开始安装过程，在xen roll安装指南中没有写这个选项，可能是开发者忽略了这一点。

（3）没有出现错误信息的话，接下来就已经开始了compute节点机的安装，这时可以在Frontend节点机上运行：

rocks-console compute-0-0-0

来监控compute节点机的安装过程，安装完以后会自动重启compute-0-0-0节点机。

（4）重复上面的方法，可以安装多个虚拟节点机。

（5）等所有虚拟节点机安装并启动以后，可以在Frontend节点机上运行：

virt-manager

来监控和管理各个虚拟机的运行情况。

至此，安装就大功告成了，至于其他的Rocks集群的其他管理命令和工具可以参照其官方网站的用户手册：

http://www.rocksclusters.org/wordpress/?page_id=4



另外，还有几个在管理Rocks集群时需要注意的地方：

1、在每次运行完rocks的一些命令修改了数据库配置信息后，比如删除compute节点机，都要再运行：

rocks sync config

来将更新后的数据库信息写入到节点机的系统配置文件中，否则在运行其他管理命令时会遇到一些莫名的错误。

2、在已有Frontend节点机上安装service-pack或添加其他新的roll的方法：

（1） wget http://www.rocksclusters.org/ftp-site/pub/rocks/rocks-5.2/linux/service-pack-5.2.2-1.i386.disk1.iso

其中wget后面的部分为在线安装时的roll文件。

<>
< init # finish-install-sp.sh tmp sh bash | service-pack roll run rocks distro create install export cd version="5.2.2" enable disable service-pack*iso add>
rocks run host '/boot/kickstart/cluster-kickstart'

2、如果集群中的节点机需要重新安装，可以在这个节点机上运行：


/boot/kickstart/cluster-kickstart


来重装系统。或者可以在Frontend节点机上运行：


rocks run host '/boot/kickstart/cluster-kickstart'


来重新安装所有的compute节点机。


3、如果想重装集群中所有的compute节点机，并在重装完以后让这些节点机继续执行由于重装而中断的计算任务，可以通过SGE控制来实现，运行：


/opt/gridengine/examples/jobs/sge-reinstall.sh


4、安装完Frontend节点机以后，centos不能正确显示中文字体，这时可以运行：


yum install fonts-chinese


来在线安装中文字体，当然这时机器必须是与互联网连通的，安装完以后将系统语言修改为中文就可以了。

cufflink ucsc genome

It works using ucsc genome gft file

1) cufflinks
/home/xusheng/Biosoft/cufflinks-1.0.3.Linux_x86_64/cufflinks -p 16 -g /home/xusheng/GenomeSequence/Human/refGene.gtf -o AD1 AD1_accepted_hits.bam

2) cuffcompare
/home/xusheng/Biosoft/cufflinks-1.0.3.Linux_x86_64/cuffcompare -V -r /home/xusheng/GenomeSequence/Human/refGene.gtf -s /home/xusheng/GenomeSequence/Human/hg19.fa -R -o s7_s8_refGene AD1/transcripts.gtf Nor3/transcripts.gtf

3. cuffdiff

home/xusheng/Biosoft/cufflinks-1.0.3.Linux_x86_64/cuffdiff -p 16 s7_s8_refGene.combined.gtf AD1/AD1_accepted_hits.bam Nor3/Nor3_accepted_hits.bam

Cufflinks problem

After running tophat, we get accepted_hits.bam file

1. /home/xusheng/Biosoft/cufflinks-1.0.3.Linux_x86_64/cufflinks -p 16 -g /home/xusheng/GenomeSequence/Human/Homo_sapiens.GRCh37.63.gtf Nor3_accepted_hits.bam

/home/xusheng/Biosoft/cufflinks-1.0.3.Linux_x86_64/cufflinks -p 16 -g Genes_July_2010_hg19.gff AD3_accepted_hits.bam (gff file only contain gene info)

remember: 1. use the UCSC gff file or ENSEMBL gff file; 2. specify the sequence file when performing the cuffcompare analysis, as the following command; 3. use -g rather than -G


3. cuffcompare step
/home/xusheng/Biosoft/cufflinks-1.0.3.Linux_x86_64/cuffcompare -o adnor_compare -s /home/xusheng/GenomeSequence/Human/Ensembl/Hs19_Ensembl.fa -r /home/xusheng/GenomeSequence/Human/Homo_sapiens.GRCh37.63.gtf AD1/AD1_transcripts.gtf AD3/AD3_transcripts.gtf AD4/AD4_transcripts.gtf AD6/AD6_transcripts.gtf AD7/AD7_transcripts.gtf Nor3/Nor3_transcripts.gtf Nor4/Nor4_transcripts.gtf Nor5/Nor5_transcripts.gtf Nor7/Nor7_transcripts.gtf

good links:

http://dingo.ucsf.edu/twiki/bin/view/Cores/BioinformaticsCore/RnaSeqAnalysis


tophat-cufflinks-cuffcompare-cuffdiff (2011-03-29 15:28:55)转载
标签： 杂谈 分类： Nextgenerationsequencing
1. Creation of a genome index (needed for Bowtie and TopHat to operate)
bowtie-build genome.fa bowtie_output/genome
2. Alignments of RNA-seq reads to the genome with identification of exon-exon splice functions using TopHat
tophat -o tophat_output/ bowtie_output/genome reads.fastq
3. Assembly of RNA-seq reads into transcripts and report of the RPKM/FPKM using Cufflinks
cufflinks -o cufflinks_output/ --reference-seq genome.fa tophat_output/accepted_hits.bam
4. Comparison of assembled transcripts with a reference annotation, or across experiments, using Cuffcompare
cuffcompare -o cuffcompare_output -r genome.gtf -R cufflinks_output/transcripts.gtf
5. Finding significant changes in transcript expression, splicing, and promoter use using Cuffdiff
cuffdiff -o cuffdiff_output/ cufflinks_output/transcripts.gtf tophat_output/accepted_hits1a.bam,tophat_output/accepted_hits1b.bam tophat_output/accepted_hits2.bam
(in this example, accepted_hits1a.bam and accepted_hits1b.sam are replicates, while accepted_hits2.bam is another experimental condition)

IP mac

IP: 128.169.4.219

在vi中如何删除^M

在vi中如何删除^M
2011-06-21 12:34:00 发表于互联网事 本文链接: 在vi中如何删除^M
如果你是用某些windows编辑器编辑的文件，很可能就会产生多余的^M转义符合，当然这些符号在wxin下是看不到的。

但是某天，你在linux下，用vim打开，发现文件中很多^M，那是多么的恶心啊，如下图：



这个时候就要vim的替换功能了，替换命令很简单 :%s/^M//g 也就是将全部的^M替换成空。

但是你可能会发现，如果你直接用键盘输入^M，是替换不成功滴，所以这里有个小技巧，^M不是直接输出来的，而是要用 ctrl+v 接下来ctrl+m 才能出现真正的文件中的^M.

上面这行，才是这篇日志记录的主要目的，因为我每次想替代^M总会忘记如何打出^M.

黄艺博在2009年武汉市总队长竞选活动上的讲话

黄艺博在2009年武汉市总队长竞选活动上的讲话

　　同志们：　

　　今天，我们在这里召开的2009年武汉市总队长竞选活动，我认为，这个会议是十分必要的，召开的相当及时。通过这个会议的召开，对于少先队工作的开展，具有十分重要的指导意义。对于刚才徐鎏杨同学的重要讲话，以及众多同学踊跃的发言，我认为，讲的非常好，非常深刻。希望在座的同志，认真领会，深刻理解。主要还是徐鎏杨领导的重要讲话，希望各位回去后，要及时向同志们传达，认真领会徐鎏杨领导的重要讲话精神，抓好落实、真抓实干，推动少先队工作的顺利开展，努力开创少先队工作的新局面。　　在这里，对于少先队的工作，我想提几点补充意见：

　　一、对于少先队的工作，我们要从思想上提高认识，充分领会少先队工作的重要性和必要性。

　　目前，少先队工作已经开创了很大的局面，获得了很大的成绩，这是有目共睹的吗！但是，还是要从深度和广度上更加推进少先队工作。我看，最重要的一点是：提高认识！是的，提高认识。各级领导要充分领会少先队工作的重要性和必要性，各级党组织要加强关于少先队工作的宣传力度，形成上下“齐抓共管”的局面，只有这样，少先队工作才能更上层楼。

　　二、 对于少先队工作，要加强落实，要把工作落到实处。

　　目前，有个别同志、个别部门，存在一个很不好的现象，就是：热衷于搞形式主义，热衷于开大会，传达文件。

　　当然，开大会嘛，还是有必要的，上传下达也是必须的。但是，光是讲空话、打官腔，是远远不够的。对少先队的工作，我们要真抓实干，加强落实。各级领导要把深入基层的工作，列入日常议事日程，要具体部署，认真执行。各级领导要为工作，创造必要的物质条件和舆论环境，扎扎实实推动工作的开展。要抓出实效，抓出成绩。

　　三、 要加强联系，搞好协调工作。

　　历史证明：钢铁一样的团结，是我武汉少先队消除一切困难的有力武器。关于工作也一样，各级领导要加强协调工作，要把上下、左右、各方面、各环节有机结合起来，步调一致地推进少先队工作的开展。目前，有些部门，遇事推诿、互相扯皮，这种官僚主义的作风是十分要不得的！这种作风，轻则导致工作效率降低，重则影响党和政府的威信。我们要坚决铲除这种官僚主义作风。有关职能部门的同志，在开展深入基层工作中，要加强联系，搞好协调工作，务必把工作做好。

　　四、要在实践中探索少先队工作与市场经济有机结合的新路子。

　　少先队工作与市场经济有没有关系，我看是大有关系。市场经济是一场深刻的社会变革，它的影响将波及社会生活的每一个领域，少先队的工作也不例外，它必然会受市场经济的影响，因此，如何适应市场经济的要求，如何和市场经济有机结合起来，希望大家认真的思考一下，去探索一下，这是十分有意义的。

　　五、参与深入基层工作的同志，要有自豪感和责任感。

　　同志们，对于少先队的工作，政府是非常重视的；各级组织也投入大量的人力、物力、财力来推进这项工作的开展。同志们，你们承担的是中国的未来，是肩负了各级组织对你们的殷切希望的，希望你们要脚踏实地、同心同德、努力工作，在各自的岗位上为社会主义建设，为改革开放，添砖加瓦！

　　以上五点意见，供同志们参考。总之，希望大家要振奋精神，多干实事，少说空话，开拓进取，努力开创少先队工作的新局面。

　　我的讲话完了。最后祝同志们五一劳动节愉快，工作顺利。为中华之崛起而奋斗！

　　谢谢各位！

Usage examples from the wild:


============================================================
Example 1
Summary: Compute the "spanning" or "physical" coverage of my paired-end
dataset on the reference genome. I will use hg18.
Contributor: Aaron Quinlan
Date: 06-July-2010
Tools used: samtools, bamToBed, genomeCoverageBed
============================================================
Workflow:
# 1. Sort the BAM files by read ID so that step 2 can correctly compute the "span" of the
# paired-end fragment.
samtools sort -n -m 50000000000 sample.sorted.dup.bam sample.namesorted.dup

# 2. Create BED from the spanning cooridnates and sort them by chrom for genomeCoverageBed.
# NOTE: Using -bedpe format from bamToBed so that I have the starts/ends for each mate.
# NOTE: Only computing coverage based upon "properly-paired" reads. (i.e., -f 0x2)
samtools view -uf 0x2 sample.namesorted.dup.bam | bamToBed -i stdin -bedpe | awk '{OFS="\t"; print $1,$2,$6,$7,$8}' | sort -k 1,1 > sample.namesorted.dup.bam.bed

# 3. Compute the ___spanning___ coverage histogram for the sample
# NOTE: We are creating a histogram where we lump all positions with depth > 60 into a single bin (-max 60)
genomeCoverageBed -i sample.namesorted.dup.bam.bed -g human.hg18.genome -max 60 > sample.spanning.cov




============================================================
Example 2
Summary: Convert a BAM into a bedgraph or a bigWig.
Contributor: Davide Cittaro
Date: 08-July-2010
Tools used: bamToBed, slopBed, genomeCoverageBed, samtools, wigToBigWig
============================================================
Workflow:
# 1. Remove duplicates from the BAM file. I'm using samtools, Picard may be also used. I'm assuming it's a single read file
samtools rmdup -s sample.bam sample.nodup.bam

# 2. Get chromosome sizes. You may get this from UCSC or from BAM headers
samtools view -H sample.nodup.bam | awk '/@SQ/ {OFS="\t"; gsub("SN:", "", $2); gsub("LN:", "", $3); print $2, $3}' > chromsizes.tab

# 3a. Convert the BAM file into a gzipped bedgraph
# NOTE: You should define the extension size (ssize) as the average library size minus the read length. There's no need to be precise here, if you have library at 200 bp and 36 bp reads, you can set ssize in a range between 160-180 bp
bamTobed -i sample.nodup.bam | slopBed -i stdin -g chromsizes.tab -s -l 0 -r ${ssize} | genomeCoverageBed -g chromsizes.tab -i stdin -bg | gzip -c - > sample.bdg.gz

# 3b. Convert the BAM file into a bigWig file
bamTobed -i sample.nodup.bam | slopBed -i stdin -g chromsizes.tab -s -l 0 -r ${ssize} | genomeCoverageBed -g chromsizes.tab -i stdin -bg | wigToBigWig stdin chromsizes.tab sample.bw



============================================================
Example 3
Summary: A slightly different way (than ex. 2) to create BigWig files from BAM
Contributor: Assaf Gordon
Date: 09-July-2010
Tools used: samtools, genomeCoverageBed, bedGraphToBigWig
More Information : http://cancan.cshl.edu/labmembers/gordon/files/viz2.pdf
============================================================
Workflow:
# 1. Convert SAM to BAM
samtools view -S -b -o sample.bam sample.sam

# 2. Sort the BAM file
samtools sort sample.bam sample.sorted

# 3. Create BedGraph coverage file
genomeCoverageBed -bg -ibam sample.sorted.bam -g chromsizes.txt > sample.bedgraph

# 4. Convert the BedGraph file to BigWig
bedGraphToBigWig sample.bedgraph chromsizes.txt sample.bw


============================================================
Example 4
Summary: Create two coverage tracks - one per strand
Contributor: Assaf Gordon
Date: 09-July-2010
Tools used: samtools, genomeCoverageBed, bedGraphToBigWig, awk
More Information : http://cancan.cshl.edu/labmembers/gordon/files/viz2.pdf
============================================================
Workflow:
# 1. Convert SAM to BAM
samtools view -S -b -o sample.bam sample.sam

# 2. Sort the BAM file
samtools sort sample.bam sample.sorted

# 3. Create BedGraph coverage file - for positive strand
genomeCoverageBed -bg -ibam sample.sorted.bam -g chromsizes.txt -strand + > sample.pos.bedgraph

# 4. Create BedGraph coverage file - for negative strand
# use AWK to negate the coverage values
genomeCoverageBed -bg -ibam sample.sorted.bam -g chromsizes.txt -strand - | awk -v OFS="\t" '{ $4 = - $4 ; print $0 }' > sample.neg.bedgraph

# 5. Create BigWig tracks from both bedgraph files
bedGraphToBigWig sample.pos.bedgraph chromsizes.txt sample.pos.bw
bedGraphToBigWig sample.pos.bedgraph chromsizes.txt sample.neg.bw

# 6. In the UCSC Genome Browser tracks, add both tracks, with different color values:
track name="Positive Strand" color=0,0,255 type=bigWig bigDataUrl=http://myserver.edu/SAMPLE.POS.BW
track name="Negative Strand" color=255,0,0 type=bigWig bigDataUrl=http://myserver.edu/SAMPLE.NEG.BW






============================================================
Example 5
Summary: Create two coverage tracks - for exons only
Contributor: Assaf Gordon
Date: 09-July-2010
Tools used: samtools, genomeCoverageBed, bedGraphToBigWig
More Information : http://cancan.cshl.edu/labmembers/gordon/files/viz2.pdf
============================================================
Workflow:
# 1. Convert SAM to BAM
samtools view -S -b -o sample.bam sample.sam

# 2. Sort the BAM file
samtools sort sample.bam sample.sorted

# 3. Create BedGraph coverage file - calculating the coverage of exons only
genomeCoverageBed -bg -ibam sample.sorted.bam -g chromsizes.txt -split > sample.bedgraph

# 4. Create BigWig tracks from the bedgraph files
bedGraphToBigWig sample.bedgraph chromsizes.txt sample.bw






=============================================================
Example 6
Summary: Find all genomic regions where all samples have >=X coverage
Contributor: Aaron Quinlan
Date: 18-July-2010
Tools used: samtools, intersectBed, genomeCoverageBed, mergeBed, awk
=============================================================
Workflow:
##############################################################
# NOTE: The example below is for three samples in the interest
# of clarity/brevity. It could be extended to many more
# samples with a BASH or Perl/Python script.
##############################################################

# 1. Let's assume you have three human samples, 1, 2, and 3.
# Let's also assume that you have aligned the sequences for
# these samples and they are now in BAM format.
ls
1.bam 2.bam 3.bam

# 1a. These BAM files contain both concordant and discordant reads.
# Also, they contain alignments with _all_ mapping qualities.
# Let's cull the files to solely those alignments that are
# "properly-paired" (0x2) and have a mapping quality >=10.
samtools view -q 10 -f 0x2 -b 1.bam > 1.clean.bam
samtools view -q 10 -f 0x2 -b 2.bam > 2.clean.bam
samtools view -q 10 -f 0x2 -b 3.bam > 3.clean.bam

# 2. First create a BEDGRAPH of the coverage for each sample.
genomeCoverageBed -ibam 1.clean.bam -g human.hg19.genome -bga > 1.clean.bam.bedg
genomeCoverageBed -ibam 2.clean.bam -g human.hg19.genome -bga > 2.clean.bam.bedg
genomeCoverageBed -ibam 3.clean.bam -g human.hg19.genome -bga > 3.clean.bam.bedg

# 3. Now merge all intervals with sufficient depth (say >=6X).
# This will create a BED file.
awk '$4 >= 6' 1.clean.bam.bedg | mergeBed -i stdin > 1.clean.bam.gte6.bed
awk '$4 >= 6' 2.clean.bam.bedg | mergeBed -i stdin > 2.clean.bam.gte6.bed
awk '$4 >= 6' 3.clean.bam.bedg | mergeBed -i stdin > 3.clean.bam.gte6.bed

# 4. Intersect them all to get the set of intervals having depth >= 6 in _all_ 3.
intersectBed -a 1.clean.bam.gte6.bed -b 2.clean.bam.gte6.bed | \
intersectBed -a stdin -b 3.clean.bam.gte6.bed > all.clean.bam.gte6.bed

# 5. What is the total number of base pairs having >=6X in _all_ 3 samples?
awk '{sum += $3-$2} END{print sum}' all.clean.bam.gte6.bed

# 6. One could now filter variant calls against this BED file to restrict
# calls to regions that you believe to be of sufficient quality, etc.