小宝的小天地吧 关注:77贴子:7,906


1楼2014-07-30 16:19回复
    什么东东?


    来自Android客户端2楼2014-08-02 09:20
    收起回复


      3楼2014-08-05 18:58
      回复


        4楼2014-08-26 13:04
        回复


          9楼2014-08-29 14:48
          回复
            perl统计文件行数


            10楼2014-08-29 14:48
            回复
              #!/usr/bin/perl -w
              system("/usr/local/bin/bowtie -p 10 -a -v 0 -f /home/databases/maize/v3 EcoR1-6.fa bowtie_out.txt ");
              my $count=0;
              open(IN,"bowtie_out.txt");
              while (<IN>){
              $count++;
              }
              my $ecor1_num=$count/2;
              close IN;
              print "EcoR1 locus number is $ecor1_num \n";


              11楼2014-08-31 20:21
              回复
                bowtie -p 10 -a -v 0 -f /home/databases/maize/v3 EcoR1.fa EcoR1_bowtie_out.txt
                #!/usr/bin/perl
                use strict;
                use warnings;
                use Getopt::Long;
                use File::Spec;
                my $small_grep_file;
                my $large_file;
                my $output;
                # Get user input
                # usage
                my $usage = "
                Usage: $0 [-ref|
                large.file] [-g|grep.file] [-o|out.file]
                Algorithm options:
                -ref large file
                -g file which is used for grepping
                -o out.file
                Note: Please read the man page for detailed
                description of the command line and options.
                \n";
                #if not input files provided, quit and print usage information
                if( !$ARGV[0] ){ die "$usage"; }
                GetOptions(
                "g|
                small_grep_file=s" => \$small_grep_file,
                "ref|large_file=s" => \$large_file,
                "o|output=s" => \$output,
                );
                open(IN,$small_grep_file) or die "$!\n";
                my $name;
                while(chomp($name=<IN>)){
                system(" cat $large_file |grep $name >>$output ");
                }
                close IN;
                exit;


                12楼2014-08-31 21:31
                回复
                  Perl不会解析单引号中的内容,但是会解析双引号中的。如果将变量放在单引号中,Perl仅仅会认为它是用户要显示字符(\’和\\除外的转义字符也不会解析),但是如果将其放在双引号的字符串里,它将被解析为一个变量。而且Perl还会解析变量字符串里的特殊字符。即使用单引号表示字符串时可以不用\作为转义字符,例如:$str = ‘This is a string’;print ‘The String is $str’;输出如下:The String is $strPerl还提供了两个函数由于引用字符串:q和qq。q函数的功能和单引号类似,qq函数的功能和双引号类似。这两个函数的主要目的是使用户在不用使用\’、\\和\”等特殊字符的情况下,就能在字符串中使用单双引号。


                  13楼2014-08-31 21:58
                  回复
                    对应11的显示结果;
                    # reads processed: 1
                    # reads with at least one reported alignment: 1 (100.00%)
                    # reads that failed to align: 0 (0.00%)
                    Reported 985460 alignments to 1 output stream(s)
                    EcoR1 locus number is 492730


                    14楼2014-08-31 22:01
                    回复
                      Bowtie2使用方法与参数详细介绍
                      Bowtie2 -q --phred33 --sensitive --end-to-end -I 0 -X 500 --fr --un unpaired --al aligned \--un-conc unconc --al-conc alconc -p 6 --reorder -x{-1-2| -U} -S []
                      用法:bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r>} -S [<hit>]
                      必须参数:-x <bt2-idx> 由bowtie2-build所生成的索引文件的前缀。首先 在当前目录搜寻,然后在环境变量 BOWTIE2_INDEXES 中制定的文件夹中搜寻。-1 <m1> 双末端测寻对应的文件1。可以为多个文件,并用逗号分开;多个文件必须和 -2 <m2> 中制定的文件一一对应。比如:"-1 flyA_1.fq,flyB_1.fq -2 flyA_2.fq,flyB_2.fq". 测序文件中的reads的长度可以不一样。-2 <m2> 双末端测寻对应的文件2.-U <r> 非双末端测寻对应的文件。可以为多个文件,并用逗号分开。测序文件中的reads的长度可以不一样。-S <hit> 所生成的SAM格式的文件前缀。默认是输入到标准输出。
                      以下是可选参数:输入参数-q 输入的文件为FASTQ格式文件,此项为默认值。-qseq 输入的文件为QSEQ格式文件。-f 输入的文件为FASTA格式文件。选择此项时,表示--ignore-quals也被选择了。-r 输入的文件中,每一行代表一条序列,没有序列名和测序质量等。选择此项时,表示--ignore-quals也被选择了。-c 后直接为比对的reads序列,而不是包含序列的文件名。序列间用逗号隔开。选择此项时,表示—ignore-quals也被选择了。-s/--skip <int> input的reads中,跳过前<int>个reads或者pairs。-u/--qupto <int> 只比对前<int>个reads或者pairs(在跳过前<int>个reads或者pairs后)。Default: no limit.-5/--trim5 <int> 剪掉5*端<int>长度的碱基,再用于比对。(default: 0).-3/--trim3 <int> 剪掉3*端<int>长度的碱基,再用于比对。(default: 0).--phred33 输入的碱基质量等于ASCII码值加上33. 在最近的illumina pipiline中得以运用。--phred64 输入的碱基质量等于ASCII码值加上64.--solexa-quals 将Solexa的碱基质量转换为Phred。在老的GA Pipeline版本中得以运用。Default: off. --int-quals 输入文件中的碱基质量为用“ ”分隔的数值,而不是ASCII码。比如 40 40 30 40...。Default: off.


                      16楼2014-08-31 22:59
                      回复
                        接上链接
                        http://www.plob.org/2012/09/02/4540.html


                        18楼2014-08-31 23:00
                        回复
                          bowtie
                          http://www.plob.org/2011/12/13/932.html


                          19楼2014-08-31 23:02
                          回复
                            http://fhqdddddd.blog.163.com/blog/static/18699154201422735644352/
                            bowtie与bowtie2 对比


                            20楼2014-08-31 23:03
                            回复
                              bowtie 比对
                              http://bowtie-bio.sourceforge.net/index.shtml
                              http://www.ncrna.net/bowtie-short-sequence-alignment-tool-detailed-solution/
                              http://www.plob.org/2011/12/13/932.html
                              常见的短序列比对工具有很多,如fasta、blast、bowtie、shrimp、soap等。每个工具都有其自身的优点,但同时也具备了一些缺点。权衡利弊,我选择bowtie作为主要的短序列比对工具。它速度很快,比对结果也容易理解。
                              现在举个例子来探讨bowtie的使用方法:现在有GENOME.fa、高通量测序数据Reads.fa,我们希望将Reads.fa比对到基因组GENOME.fa上。
                              (一)、对Reference文件(GENOME.fa)建库
                              1 bowtie-build GENOME.fa GENOME.fa
                              建库步骤可能需要1h甚至几个小时,建议在后台执行:
                              nohup bowtie-build GENOME.fa GENOME.fa &
                              (二)、将Reads.fa比对到GENOME.fa上,只能比对到正链,且匹配到基因组不多于20个不同位置,允许有1个错配
                              1 bowtie -f -a -m 20 -v 1 --al Reads_aligned --un Reads_unaligned --norc GENOME.fa Reads.fa Reads.bwt 2> log
                              注:
                              -f 指定query文件为fasta格式
                              -a 保留所有比对结果
                              -m 指定最大比对到基因组的次数
                              -v 允许最大错配数,为[0-2]
                              --al 能map到GENOME的reads,fasta格式
                              --un 不能map到GENOME的reads,fasta格式
                              --norc 不输出匹配到负链的结果;如果不想输出比对到正链的结果,则用"--nofw"。不指定该选项则正负链结果都输出
                              后面依次写上GENOME索引文件,Reads文件,输出结果文件Reads.bwt,日志文件log。
                              (三)、bowtie输出结果的说明
                              12 sample001_x75 + Chr1 12453 ATCGGCCAATTACGGACTTAA IIIIIIIIIIIIIIIIIIIII 4 9:G>T 1 2 3 4 5 6 7 8
                              1. query id
                              2. "+"表示正向match;"-"表示对query作反向互补后match
                              3. reference id
                              4. 第2列为"+"时,表示query 第一个碱基map到reference(5*->3*)上的位置,0-based(以0开始);第2列为"-"时,表示query的反向互补序列第 一 个碱基map到reference(5*->3*)上的位置,0-based(以0开始)
                              5. 如果第2列为"+",则和query序列一致;否则,和query序列反向互补
                              6. 质量文件,如果query文件为fasta格式,则无法获取质量文件,用I代替,I的数量与query序列长度一致
                              7. 当前query能map到GENOME的4个不同位置
                              8. 如果存在第8列,表示有mismatch。第8列可以分为三个部分,最左端的数字,中间的碱基为reference碱基,最右端的碱基为query碱基,下面分情况讨论:
                              第2列为"+"时:最左端的数字9表示query从5*端数起,第10个碱基为"T",而对应的reference为"G";
                              第2列为"-"时:最左端的数字9表示query先作反向互补,然后从3*端数起,第10个碱基为"T",而对应的reference为"G";


                              21楼2014-08-31 23:04
                              回复