ssh不提示将ip加入到known_hosts
需要使用ssh连接到指定的主机,然后运行一个命令后退出。但是该主机每次重新联网都会变化IP(PPOE上网).这时每次ssh的连接默认情况下会询问是否将主机ip永远加入到本地的~/.ssh/know_hosts中。
这种情况下则所需要的命令不会运行。比如运行程序:
$ ssh user@ip_address "mplayer ~/musics/audio/alarms/Beeps.mp3"
以上命令的目是主机上的声音。将此命令在远端服务器运行,来播放本地机上的一个提醒铃声。但是本地主机没有固定IP的时候,就不会有提醒。于是进行如下设置:
更改/etc/ssh/ssh_config文件
StrictHostKeyChecking no
只根据两个转录组数据分析差异表达基因
1. 数据
手头有两个转录组数据,分别是某一真菌物种(该物种没有基因组数据)的双核菌丝阶段和子实体阶段的转录组数据。测序平台是Illumina Hiseq2000;插入片段长度200bp,测序的reads长度90bp。得到的数据文件为:
/home/user/RNA-seq/mycelium_reads1.fastq /home/user/RNA-seq/mycelium_reads2.fastq /home/user/RNA-seq/fruitingbody_reads1.fastq /home/user/RNA-seq/fruitingbody_reads2.fastq
2. 数据的预处理
2.1. 去除N的比例大于5%的reads;去除低质量reads(质量值Q≤20的碱基数占整个read的50%以上);
2.2. 根据FastQC对上述过滤后的reads的质量检测,去除reads首尾各10bp的碱基,得到的预处理数据为:
/home/user/RNA-seq/clean_reads/mycelium_reads1.fastq /home/user/RNA-seq/clean_reads/mycelium_reads2.fastq /home/user/RNA-seq/clean_reads/fruitingbody_reads1.fastq /home/user/RNA-seq/clean_reads/fruitingbody_reads2.fastq
3. 使用Trinity和TGICL进行转录组的组装
3.1. 对mycelium数据和fruitingbody数据分别进行转录组的组装
$ pwd /home/user/RNA-seq/ $ mkdir assembly $ cd assembly $ Trinity.pl --jaccard_clip --seqType fq --JM 50G --SS_lib_type FR --CPU 24 --inchworm_cpu 24 --bflyCPU 24 --group_pairs_distance 500 -min_contig_length 200 --output mycelium_contig --left /home/user/RNA-seq/clean_reads/mycelium_reads1.fastq --right /home/user/RNA-seq/clean_reads/mycelium_reads2.fastq $ Trinity.pl --jaccard_clip --seqType fq --JM 50G --SS_lib_type FR --CPU 24 --inchworm_cpu 24 --bflyCPU 24 --group_pairs_distance 500 -min_contig_length 200 --output fruitingbody_contig --left /home/user/RNA-seq/clean_reads/fruitingbody_reads1.fastq --right /home/user/RNA-seq/clean_reads/fruitingbody_reads2.fastq 当然,需要将生成的两个转录组的Trinity.fasta序列按长度进行排序;统一更改的fasta 头名称;更改fasta文件名 $ cp mycelium_contig/Trinity.fasta mycelium_contigs.fasta $ cp fruitingbody_contig/Trinity.fasta fruitingbody_contigs.fasta
3.2. 使用TGICL将两个转录组序列合并
$ pwd /home/user/RNA-seq/assembly $ mkdir all_tissue_contig $ cd all_tissue_contig $ cat ../mycelium_contigs.fasta ../fruitingbody_contigs.fasta > all.contigs.fasta $ tgicl -F all.contigs.fasta tgicl生成了两个有用的文件: asm_1/contigs 和 all.contigs.fasta.single tons。其中前者是聚类后的contigs结果;后者是没有聚类的单独的contigs的序列名,需 要分别到../mycelium_contigs.fasta 和 ../fruitingbody_contigs.fasta文 件中提取出相应的序列: all.contigs.fasta.singletons.mycelium.fasta 和 all.contigs.fasta.singletons.fruitingbody.fasta $ cat asm_1/contigs all.contigs.fasta.singletons.mycelium.fasta all.contigs.fasta.singletons.fruitingbody.fasta > all_contigs.fasta 在对all_contigs.fasta进行fasta头的重命名,并将序列按长度排序;同时要得到all_ contigs.fasta中的序列和mycelium_contigs.fasta,fruitingbody_contigs. fasta中序列的对应关系。 $ cp all_contigs.fasta ../
3.3. 至此得出转录组的组装结果:
/home/user/RNA-seq/assembly/all_contigs.fasta
4. 使用cufflinks来分析差异表达基因
4.1 使用tophat来将预处理后的reads比对到转录组序列上
$ pwd /home/user/RNA-seq/ $ mkdir cufflinks $ cd cufflinks $ bowtie2-build ../all_contigs.fasta all_contigs $ tophat -o tophat_out_mycelium -r 60 --mate-std-dev 80 -p 24 --library-type fr-unstranded all.contigs /home/user/RNA-seq/clean_reads/mycelium_reads1.fastq /home/user/RNA-seq/clean_reads/mycelium_reads2.fastq $ tophat -o tophat_out_fruitingbody -r 60 --mate-std-dev 80 -p 24 --library-type fr-unstranded all.contigs /home/user/RNA-seq/clean_reads/fruitingbody_reads1.fastq /home/user/RNA-seq/clean_reads/fruitingbody_reads2.fastq
4.2 自制一个transcripts.gtf文件
由于cuffdiff运行需要一个参考序列的transcripts.gtf文件: 该文件有9列,使用tab分隔;使exon的范围为整个contig。其格式如:
AA.auricula_all_contig_1 chenlianfu exon 1 9401 . . . gene_id "A.auricula_all_contig_1"; transcript_id "A.auricula_all_contig_1"; A.auricula_all_contig_10 chenlianfu exon 1 4464 . . . gene_id "A.auricula_all_contig_10"; transcript_id "A.auricula_all_contig_10"; A.auricula_all_contig_100 chenlianfu exon 1 3090 . . . gene_id "A.auricula_all_contig_100"; transcript_id "A.auricula_all_contig_100"; A.auricula_all_contig_1000 chenlianfu exon 1 1768 . . . gene_id "A.auricula_all_contig_1000"; transcript_id "A.auricula_all_contig_1000"; A.auricula_all_contig_10000 chenlianfu exon 1 586 . . . gene_id "A.auricula_all_contig_10000"; transcript_id "A.auricula_all_contig_10000"; A.auricula_all_contig_10001 chenlianfu exon 1 586 . . . gene_id "A.auricula_all_contig_10001"; transcript_id "A.auricula_all_contig_10001";
4.3 使用cuffdiff来分析转录子的表达量和差异表达基因
$ pwd /home/user/RNA-seq/cufflinks $ cuffdiff -L mycelium,fruitingbody --library-type fr-unstranded -p 8 -o cuffdiff ./transcriptome.gtf ./tophat_out_mycelium/accepted_hits.bam ./tophat_out_fruitingbody/accepted_hits.bam
至此得出转录子的表达量数据和差异表达分析
/home/user/RNA-seq/cufflinks/cuffdiff/gene_exp.diff
基因组组装软件心得
1. SOAPdenovo 和 Velvet组装基因组是否需要对reads进行预处理?
使用soapdenovo组装基因组,需要先对reads进行预处理,去掉质量不好的reads,剪掉3’端和5’端的一些碱基;然后再进行组装。而使用velvet时候,则不需要,可能是由于velvet的cutoff比较好。
这个结论得出方法:将reads进行预处理前后分别用于组装出基因组,然后再使用bowtie2将两种reads比对到两种基因组组装结果,观察reads比对上的百分比得出。reads的预处理方法是去除低质量reads并将100bp的reads前后各减去5bp碱基。
对于soapdenovo:
未处理reads比对率 预处理reads比对率 未处理reads基因组 65.84% 66.10% 预处理reads基因组 66.68% 66.99%
对于velvet:
未处理reads比对率 预处理reads比对率 未处理reads基因组 86.41% 88.10% 预处理reads基因组 70.09% 71.62%
通过以上的对比和数据可以看出:1.velvet组装出的基因组的精确性要比soapdenovo高;2.对reads进行预处理对soapdenovo组装的基因组的精确性影响较小,预处理reads后,组装出的基因组精确性稍有提高;3.使用velvet组装基因组时,不对reads进行预处理组装出的基因组精确性高很多。
讨论:1.数据量不够多,碱基覆盖度只有20~25X, 以上结论具有一定的片面性; 2.velvet在对kmer进行cutoff的算法很棒,因而组装出的基因组精确性高。
Velvet的安装和使用
一. Velvet简介
Velvet 是 EBI 的 Daniel Zerbino and Ewan Birney 开发的利用short reads组装基因组的软件。
Velvet作为一种de novo基因组组装软件,可利用的reads包括 Illumina reads 和 454 reads。
引用文献(2008):Velvet: algorithms for de novo short read assembly using de Bruijn graphs.
二. Velvet下载与安装
make ‘CATEGORIES=57’ ‘MAXKMERLENGTH=75’ ‘BIGASSEMBLY=1’ ‘LONGSEQUENCES=1’ ‘OPENMP=1’ ‘BUNDLEDZLIB=1’
远程服务器的命令——到本地声音提醒
使用ssh连接远程服务器后。每当运行一个较长的程序时候,不想一直等待,需要该程序结束后能有声音提醒。这时候可以参考我的使用方法:
1. 下载一些提醒铃声到本地文件中
本地一个铃声的存放路径: ~/musics/audio/alarms/Beeps.mp3
2. 和远程服务器相互设置无密码的ssh登录
参考我以前的一篇文章:CentOs配置OpenSSH无密码登陆
3. 编写脚本alarm.pl放置于服务器的系统路径下
#!/usr/bin/perl use strict; my $whoami = `whoami`; $whoami =~ s/\s*//g; my $lastlog = `lastlog`; my $ip_address; if ( $lastlog =~ m/^$whoami.*?(\d+\.\d+\.\d+\.\d+)/m ) { $ip_address = $1; } else { die; } print "The ip_address is $ip_address\n"; my $commond = "ssh chenlianfu\@$ip_address \"mplayer -loop 0 ~/musics/audio/alarms/Beeps.mp3\""; print "$commond\n"; `$commond`;
4. 每次执行运算时间较长的程序时则紧跟着运行 alarm.pl 则会在程序运行完毕后在本地计算机上提醒。
linux下ISO文件的制作,刻录和挂载
1. ISO文件的刻录
手头上有iso文件,比如CentOS-6.4-x86_64-bin-DVD1.iso,需要刻录成光碟,然后用于装系统用。使用 cdrecord 命令。
$ cdrecord dev=/dev/cdrw CentOS-6.4-x86_64-bin-DVD1.iso 常用参数: -v | -verbose speed=<int> 指定刻录机速度,一般不用,因为cdrecord会自动检测最佳刻录速度 dev=<target> 光驱的位置
2. ISO文件的制作
制作ISO文件可以是用 dd 或 mkisofs 命令。
2.1 dd用法
当直接从光驱中压制ISO文件时:
$ dd if=/dev/cdrw of=output.iso bs=1024
2.2 mkisofs用法
当将一个文件或目录压制到ISO文件时:
$ mkisofs -r -o output.iso folder_name/
3. ISO文件的挂载
$ sudo mount -o loop file.iso /mnt
wordpress插件
1. Akismet
用途:博客垃圾评论 / trackback 屏蔽解决方案。效果非常好。
需要的Akismet API 密钥:c28ea9690585
2. WP-Chinese-Excerpt
用途:在主页上只显示中文摘要。开启该插件则会有效果。编辑文章的时候最好写个摘要。其原版插件名为:WP-UTF8-Excerpt
Tophat的安装与使用
一. Tophat简介
Tophat使用RNA-seq的reads数据来寻找基因的剪切点(splice junction)。该软件调用Bowtie,或Bowtie2来将reads比对到参考基因组上,分析比对结果,从而寻找出外显子之间的结合位点。
二. Tophat安装
直接下载适合于Linux x86_64的二进制文件,解压缩即可使用。
$ wget http://tophat.cbcb.umd.edu/downloads/tophat-2.0.8b.Linux_x86_64.tar.gz $ tar zxf tophat-2.0.8b.Linux_x86_64.tar.gz
前提条件当然要安装Bowtie, Bowtie2, SAM tools, Boost C++ libraries等。
三. Tophat的使用参数
使用Tohat时,bowtie2(或bowtie,下同), bowtie2-align, bowtie2-inspect, bowtie2-build 和 samtools 必须要在系统路径中。
1. 用法
$ tophat [options]* <index_base> <reads1_1[,...,readsN_1]> <reads1_2[,...,readsN_2]> 可以看出,tophat必须要的条件是比对的index数据库,以及要比对的reads。可以为多个 paired-end reads数据以逗号分开。
值得注意的:Tophat能比对的最大reads为1024bp; 能比对paired-end reads; 不能将多种不同类型的reads混合起来进行比对,这样会给出不好的结果。
如果有多种不同类型的reads进行比对,则可以:
首先,对一种类型的reads使用合适的参数运行tophat; 接着,使用bed_to_juncs将前一次的运行结果junctions.bed转换成下一次运行tophat 所的-j参数所需的junction文件; 最后,再一次使用-j参数运行tophat。
2. 常用一般参数
-h | --help -v | --version -N | --read-mismatches default: 2 丢弃不匹配碱基数超过该数目的比对结果 --read-gap-length default: 2 丢弃gap总长度超过该数目的比对结果 --read-edit-dist default: 2 丢弃read的edit distance大于该值的比对结果 --read-realign-edit-dist <int> default: "read-edit-dist" + 1 一些跨越多个exons的reads可能会被错误地比对到geneome上。Tophat有多个比对 步骤,每个比对步骤过后,比对结果中包含了edit distance的值。该参数能让Tophat对 那些edit distance的值 >= 该参数的reads重新进行比对。若设置该参数值为0,则每个 read在多个比对步骤中每次都要进行比对。这样会加大地增加比对精确性和运行时间。默认下 该参数比上一个参数的值大,则表示对reads进行重新比对。 --bowtie1 default: bowtie2 使用Bowtie1来代替Bowtie2进行比对。特别是使用colorspace reads时,因为只 有Bowtie1支持,而Bowtie2不支持。 -o | --output <string> default: ./tophat_out 输出的文件夹路径 -r | --mate-inner-dist <int> default: 50 成对的reads之间的平均inner距离。例如:fragments长度300bp,reads长度50bp , 则其inner距离为200bp,该值该设为200。 --mate-std-dev <int> default:20 inner距离的标准偏差。 -a | --min-anchor-length <int> default: 8 read的锚定长度:该参数能设定的最小值为3;锚定在junction两边的reads长度只 有都大于此值,才能用于junction的验证。 -m | --splice-mismatches <int> default: 0 对于一个剪切比对,其在锚定区能出现的最大的不匹配碱基数。 -i | --min-intron-length <int> default:70 最小的intron长度。Tophat会忽略比该长度要小的donor/acceptor pairs,认 为该区属于exon。 --I | --max-intron-length <int> default:500000 最大的intron长度。Tophat会忽略长度大于该值的donor/acceptor pairs,除 非有long read支持。 --max-insertion-length <int> defautl: 3 最大的插入长度 --max-deletion-length <int> default: 3 最大的缺失长度 --solexa-quals fastq文件使用Solexa的碱基质量格式 --solexa1.3-quals | --phred64-quals 使用Illumina GA pipeline version 1.3的碱基质量格式,即Phred64. -Q | --quals 说明是使用单独的碱基质量文件 --inter-quals 有空格隔开的整数值来代表碱基质量。当使用 -C 参数时,该参数为默认参数。 -C | --color Colorspace reads。使用这一种reads的时候命令如下: $ tophat --color --quals --bowtie1 [other options]* <colorspac e_index_base> <reads1_1[,...,readsN_1]> <reads2_1[,...,readN_2]> <quals1_1[,...,qualsN_1]> <quals1_2[,...,qualsN_2]> -p | --num-threads <int> default: 1 比对reads的线程数 -g | --max-multihits <int> default: 20 对于一个reads,可能会有多个比对结果,但tophat根据比对得分,最多保留的比对结 果数目。如果没有 --report-secondary-alignments 参数,则只会报告出最佳的比对 结果。若最佳比对结果数目超过该参数值,则只随机报告出该数目的最佳比对结果;若有 -- report-secondary-alignments 参数,则按得分顺序报告出比对结果,直至达到默认 的数目为止。 --report-secondary-alignments 是否报告additional or secondary alignments(基于比对分值AS来确定的)。 --no-discordant 对于paired reads,仅仅报告concordant mappings。 --no-mixed 对于paired reads,只报告concordant mappings 和 discordant mappi ngs。默认上,是所有的比对结果都报告。 --no-coverage-search 取消以覆盖度为基础来搜寻junctions,和下一个参数对立,该参数为默认参数。 --coverage-search 确定以覆盖度为基础来搜寻junctions。该参数能增大敏感性。 --microexon-search 使用该参数,pipeline会尝试寻找micro-exons。仅仅在reads长度>=50bp时有效。 --library-type Tophat处理的reads具有链特异性。比对结果中将会有个XS标签。一般Illumina数 据的library-type为 fr-unstranded。
3. 高级参数
--keep-tmp 保留中间文件和临时文件,对于debug有用 --keep-fasta-order 对比对结果按基因组fasta文件进行排序。该参数会使输出的SAM/BAM文件和tophat的 1.41或以前版本不兼容 --no-sort-bam 输出的BAM文件不是coordinate-sorted. --no-convert-bam 不要转换成bam格式。输出结果为sam格式。 -R | --resume <string> 从最末尾的成功完成点处,接着运行Tophat。使用方法为: $ tophat -R tophat_out -z | --zpacker default:gzip 用来对临时文件进行压缩的的压缩程序
4. Bowtie2的特别参数
使用tophat2的时候,其中的一些参数传递为bowtie2的参数,这些参数都以’b2’开头。其实,这些参数使用默认的即可。
end-to-end模式(Tophat2不能使用local alignment): --b2-very-fast --b2-fast --b2-sensitive --b2-very-sensitive 比对参数: --b2-N default: 0 --b2-L default: 20 --b2-i default: S,1,1.25 --b2-n-ceil default: L,0,0.15 --b2-gbar defaut: 4 得分参数: --b2-mp default: 6,2 --b2-np default: 1 --b2-rdg default: 5,3 --b2-rfg default: 5,3 --b2-score-min default: L,-0.6,-0.6 Effort参数 --b2-D default: 15 --b2-R default: 2
5. 融合转录子mapping
如果设定 –fusion-search 参数,则有些reads能比对到潜在的融合转录子(fusion transcripts)上。额外融合信息保存在 fusions.out 中。
--fusion-search 开启融合转录子的比对 --fusion-anchor-length default: 20 read比对到融合子的两边,每以边至少匹配的碱基数。
6. 提供的转录子的结构注释数据
值得注意的提供的GTF文中中的染色体名称和Bowtie index中的一致。这些名称是区分大小写的。
-j | --raw-juncs <.juncs file> 提供junctions文件。该文件可以使用tophat同一目录下的程序bed_to_juncs程序 来处理tophat的结果文件junctions.bed生成。 $ bed_to_juncs <junctions.bed> junctions文件是tab分隔的文件,内容为: <chrom> <leftgt> <rigth> <+/-> 其中left和right数值是0-based的junction两端的值。 --no-novel-juncs 只搜寻和GFF或junctions文件中提供的junctions想匹配的reads。如果没有 -G 或 -j 参数,则该参数无效。 -G | --GTF <GTF/GFF3 file> 提供基因模型的注释文件,GTF 2.2 或者 GFF 3 格式的文件。如果设置了该参数,Tophat 则先提取出转录子序列,然后使用Bowtie2将reads比对到提取的转录组中;只有不能比对上 的reads再比对到genome;比对上的reads再打断转变成genomic mappings;再融合新 的mappings和junctions作为最后的输出。 值得注意的是GTF/GFF文件代表chromosome和contig的第一列要和bowtie index中的 参考序列名一致。 `$ bowtie-inspect --names your_index` 命令可以获得bowt ie的index。 --transcriptome-index <dir/prefix> 是使用了 -G 参数后,Tophat提取转录子序列,然后使用bowtie2-build来建立index, 这个过程会消耗不少时间。于是,使用该参数,会将index文件生成到指定文件夹。则后续的 运用同样的index则不再需要额外耗时了。
7. 提供insertions/deletions
以下参数是使用RNA-seq数据来验证indels。
--insertions | --deletions <.juncs file> juncs文件例子: chr1 20564 20567 0-based数值。表示有20565和20566这2个碱基缺失 chr1 17491 17491 CA 表示在17491处插入了2个碱基CA --no-novel-indels 仅仅只搜寻在已给的位点的reads。
四. Tophat的输出结果
主要的结果文件是:
1. accepted_hits.bam
2. junctions.bed UCSC BED格式
3. insertions.bed 和 deletions.bed
五. 思考题
对一物种的两个样本A和B使用Illumina Hiseq2000分别进行了转录组测序,得到了结果文件A.reads1_1.fastq, A.reads1_2.fastq, B.reads1_1.fastq 和 B.reads1_2fastq。测序文库的插入片段长度为200bp,reads长度为90bp。物种的基因组文件为species.fasta。请用Tophat分析该物种的转录结合位点,indel信息?
$ bowtie2-build species.fastq species $ tophat \ --read-realign-edit-dist 0 \ -o ./tophat_out \ -r 20 \ --mate-std-dev 20 \ --coverage-search \ --microexon-search \ -p 24 \ --library-type fr-unstranded \ species \ A.reads1_1.fastq,B.reads1_1.fastq A.reads1_2.fastq B.reads1_2fastq
perl笔记
1. select函数
select FileHandl
该函数用来设置输出的文件句柄为默认文件句柄,并返回以前默认的文件句柄。最原始的默认文件句柄是 ‘main::STDOUT’ ,即标准输出。为了输出的方便,可以使用该函数指定默认的输出。
使用例子:
open FILE, '>', "job"; $oldHandle = select ( FILE ); print "This is sent to $oldHandle\n";print "This is sent to $oldHandle\n"; $origin = select ( $oldHandle ); print "This is sent to $origin\n";
2. 查找perl模块或卸载
使用以下命令列出所有perl模块的位置。eg:通过grep GD,来找出GD模块所处位置,然后直接rm掉即可。
$ perl -MFile::Find -le 'finddepth({wanted=>sub{print $_ if/\.pm$/},no_chdir=>1},@INC)' | grep 'GD'
/usr/local/lib64/perl5/GD.pm
/usr/local/lib64/perl5/GD/Simple.pm
/usr/local/lib64/perl5/GD/Polygon.pm
/usr/local/lib64/perl5/GD/Image.pm
/usr/local/lib64/perl5/GD/Group.pm
/usr/local/lib64/perl5/GD/Polyline.pm
/usr/local/share/perl5/GD/SVG.pm
/usr/local/share/perl5/Bio/Graphics/GDWrapper.pm
/usr/lib64/perl5/GDBM_File.pm
3. 十进制转二进制
$a = sprintf("%b",$a)+0;
$ perl -MFile::Find -le 'finddepth({wanted=>sub{print $_ if/\.pm$/},no_chdir=>1},@INC)' | grep 'GD' /usr/local/lib64/perl5/GD.pm /usr/local/lib64/perl5/GD/Simple.pm /usr/local/lib64/perl5/GD/Polygon.pm /usr/local/lib64/perl5/GD/Image.pm /usr/local/lib64/perl5/GD/Group.pm /usr/local/lib64/perl5/GD/Polyline.pm /usr/local/share/perl5/GD/SVG.pm /usr/local/share/perl5/Bio/Graphics/GDWrapper.pm /usr/lib64/perl5/GDBM_File.pm
3. 十进制转二进制
$a = sprintf("%b",$a)+0;