文章目录
- 一、在任意文件夹下面创建形如 1/2/3/4/5/6/7/8/9 格式的文件夹系列
- 二、在创建好的文件夹下面,比如我的是 /Users/jimmy/tmp/1/2/3/4/5/6/7/8/9 ,里面创建文本文件 me.txt
- 三、在文本文件 me.txt 里面输入内容:
- 四 、 删除上面创建的文件夹 1/2/3/4/5/6/7/8/9 及文本文件 me.txt
- 五、在任意文件夹下面创建 folder1~5这5个文件夹,然后每个文件夹下面继续创建 folder1~5这5个文件夹,效果如下:
- 六、在第五题创建的每一个文件夹下面都 创建第二题文本文件 me.txt ,内容也要一样。
- 七,再次删除掉前面几个步骤建立的文件夹及文件
- 八、下载 http://www.biotrainee.com/jmzeng/igv/test.bed 文件,后在里面选择含有 H3K4me3 的那一行是第几行,该文件总共有几行。
- 九、下载 http://www.biotrainee.com/jmzeng/rmDuplicate.zip 文件,并且解压,查看里面的文件夹结构
- 十、打开第九题解压的文件,进入 rmDuplicate/samtools/single 文件夹里面,查看后缀为 .sam 的文件,搞清楚 生物信息学里面的SAM/BAM 定义是什么。
- 十一、安装 samtools 软件
- 十二、打开 后缀为BAM 的文件,找到产生该文件的命令。 提示一下命令是:
- 十三题、根据上面的命令,找到我使用的参考基因组 /home/jianmingzeng/reference/index/bowtie/hg38 具体有多少条染色体
- 十五题、重新打开 rmDuplicate/samtools/paired 文件夹下面的后缀为BAM 的文件,再次查看第二列,并且统计
- 十六题、下载 http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip 文件,并且解压,查看里面的文件夹结构, 这个文件有2.3M,注意留心下载时间及下载速度。
- 十七题、解压 sickle-results/single_tmp_fastqc.zip 文件,并且进入解压后的文件夹,找到 fastqc_data.txt 文件,并且搜索该文本文件以 >>开头的有多少行?
- 十八题、下载 http://www.biotrainee.com/jmzeng/tmp/hg38.tss 文件,去NCBI找到TP53/BRCA1等自己感兴趣的基因对应的 refseq数据库 ID,然后找到它们的hg38.tss 文件的哪一行。
- 十九题、解析hg38.tss 文件,统计每条染色体的基因个数。
- 二十题、解析hg38.tss 文件,统计NM和NR开头的序列,了解NM和NR开头的含义。
一、在任意文件夹下面创建形如 1/2/3/4/5/6/7/8/9 格式的文件夹系列
[sunchengquan 21:42:32 ~/test]
$ mkdir -p 1/2/3/4/5/6/7/8/9
[sunchengquan 21:42:55 ~/test]
$ ls
000.csv chrY_SNP_pos_id.vcf id rs.csv sun
1 encode.py info rs_split_v1.txt test.bed
123.csv gff merge_file.sh sk.sh tmp
allid grasp_ncbi_snp01.py name.txt socket_client.py tmp.sh
arg1.txt grasp_ncbi_snp.py output_file socket_server.py
arg.txt hum_name.txt result.csv sum_name.csv
[sunchengquan 21:43:04 ~/test]
$ cd 1/2/3/4/5/6/7/8/9
[sunchengquan 21:43:46 ~/test/1/2/3/4/5/6/7/8/9]
$ pwd
/home/sunchengquan/test/1/2/3/4/5/6/7/8/9
二、在创建好的文件夹下面,比如我的是 /Users/jimmy/tmp/1/2/3/4/5/6/7/8/9 ,里面创建文本文件 me.txt
[sunchengquan 21:43:50 ~/test/1/2/3/4/5/6/7/8/9]
$ touch me.txt
[sunchengquan 21:44:51 ~/test/1/2/3/4/5/6/7/8/9]
$ ls
me.txt
三、在文本文件 me.txt 里面输入内容:
Go to: http://www.biotrainee.com/
I love bioinfomatics.
And you ?
[sunchengquan 15:29:52 ~/test]
$ cat > me.txt
Go to: http://www.biotrainee.com/
I love bioinfomatics.
And you ?
^C
[sunchengquan 15:36:18 ~/test]
$ cat me.txt
Go to: http://www.biotrainee.com/
I love bioinfomatics.
And you ?
[sunchengquan 21:44:53 ~/test/1/2/3/4/5/6/7/8/9]
$ vi me.txt
[sunchengquan 21:45:56 ~/test/1/2/3/4/5/6/7/8/9]
$ cat me.txt
Go to: http://www.biotrainee.com/
I love bioinfomatics.
And you ?
四 、 删除上面创建的文件夹 1/2/3/4/5/6/7/8/9 及文本文件 me.txt
[sunchengquan 21:46:01 ~/test/1/2/3/4/5/6/7/8/9]
$ rm -r ~/test/1
五、在任意文件夹下面创建 folder1~5这5个文件夹,然后每个文件夹下面继续创建 folder1~5这5个文件夹,效果如下:
[sunchengquan 15:46:04 ~/test/sun]
$ mkdir -p folder_{1..5}/folder_{1..5}
[sunchengquan 15:46:10 ~/test/sun]
$ ls *
folder_1:
folder_1 folder_2 folder_3 folder_4 folder_5
folder_2:
folder_1 folder_2 folder_3 folder_4 folder_5
folder_3:
folder_1 folder_2 folder_3 folder_4 folder_5
folder_4:
folder_1 folder_2 folder_3 folder_4 folder_5
folder_5:
folder_1 folder_2 folder_3 folder_4 folder_5
六、在第五题创建的每一个文件夹下面都 创建第二题文本文件 me.txt ,内容也要一样。
[sunchengquan 15:46:24 ~/test/sun]
$ touch folder_{1..5}/folder_{1..5}/me.txt
简洁明了的方法
[sunchengquan 15:56:34 ~/test/sun]
$ cat > me.txt
sun
cheng
quan
^C
[sunchengquan 17:31:34 ~/test/sun]
$ echo folder_{1..5}/folder_{1..5} |xargs -n1 cp -v me.txt
"me.txt" -> "folder_1/folder_1/me.txt"
"me.txt" -> "folder_1/folder_2/me.txt"
"me.txt" -> "folder_1/folder_3/me.txt"
"me.txt" -> "folder_1/folder_4/me.txt"
"me.txt" -> "folder_1/folder_5/me.txt"
"me.txt" -> "folder_2/folder_1/me.txt"
"me.txt" -> "folder_2/folder_2/me.txt"
"me.txt" -> "folder_2/folder_3/me.txt"
"me.txt" -> "folder_2/folder_4/me.txt"
"me.txt" -> "folder_2/folder_5/me.txt"
"me.txt" -> "folder_3/folder_1/me.txt"
"me.txt" -> "folder_3/folder_2/me.txt"
"me.txt" -> "folder_3/folder_3/me.txt"
"me.txt" -> "folder_3/folder_4/me.txt"
"me.txt" -> "folder_3/folder_5/me.txt"
"me.txt" -> "folder_4/folder_1/me.txt"
"me.txt" -> "folder_4/folder_2/me.txt"
"me.txt" -> "folder_4/folder_3/me.txt"
"me.txt" -> "folder_4/folder_4/me.txt"
"me.txt" -> "folder_4/folder_5/me.txt"
"me.txt" -> "folder_5/folder_1/me.txt"
"me.txt" -> "folder_5/folder_2/me.txt"
"me.txt" -> "folder_5/folder_3/me.txt"
"me.txt" -> "folder_5/folder_4/me.txt"
"me.txt" -> "folder_5/folder_5/me.txt”
[sunchengquan 17:33:29 ~/test/sun]
$ cat folder_5/folder_5/me.txt
sun
cheng
quan
麻烦琐碎的方法
#!/usr/bin/bash
mkdir -p ~/test/sun/folder_{1..5}/folder_{1..5}
for i in {1..5};do
cd ~/test/sun/folder_$i
for a in {1..5};do
cd ~/test/sun/folder_$i/folder_$a
echo -e "sun\ncheng\nquan\n I love you" >me.txt
done
done
七,再次删除掉前面几个步骤建立的文件夹及文件
[sunchengquan 17:36:39 ~/test/sun]
$ rm -r folder*
[sunchengquan 17:36:47 ~/test/sun]
$ ll
总用量 0
八、下载 http://www.biotrainee.com/jmzeng/igv/test.bed 文件,后在里面选择含有 H3K4me3 的那一行是第几行,该文件总共有几行。
[sunchengquan 17:36:49 ~/test/sun]
$ wget -c http://www.biotrainee.com/jmzeng/igv/test.bed
[sunchengquan 17:39:29 ~/test/sun]
$ ll
总用量 4.0K
-rw-r----- 1 sunchengquan sunchengquan 3.1K 5月 18 2017 test.bed
[sunchengquan 17:41:21 ~/test/sun]
$ grep -n 'H3K4me3' test.bed
8:chr1 9810 10438 ID=SRX387603;Name=H3K4me3%20(@%20HMLE);Title=GSM1280527:%20HMLE%20Twist3D%20H3K4me3%20rep2%3B%20Homo%20sapiens%3B%20ChIP-Seq;Cell%20group=Breast;
source_name=HMLE_Twist3D_H3K4me3;cell%20type=human%20mammary%20epithelial%20cells;transfected%20with=Twist1;culture%20type=sphere;chip%20antibody=H3K4me3;chip%20antibody%20vendor=Millipore; 222 9810 10438 0,226,255
[sunchengquan 17:42:59 ~/test/sun]
$ wc -l test.bed
10 test.bed
九、下载 http://www.biotrainee.com/jmzeng/rmDuplicate.zip 文件,并且解压,查看里面的文件夹结构
[sunchengquan 17:44:31 ~/test/sun]
$ wget -c http://www.biotrainee.com/jmzeng/rmDuplicate.zip
[sunchengquan 17:45:33 ~/test/sun]
$ ll
总用量 112K
-rw-r----- 1 sunchengquan sunchengquan 103K 11月 12 2016 rmDuplicate.zip
-rw-r----- 1 sunchengquan sunchengquan 3.1K 5月 18 2017 test.bed
[sunchengquan 17:45:48 ~/test/sun]
$ unzip rmDuplicate.zip
[sunchengquan 17:45:59 ~/test/sun]
$ ll
总用量 116K
drwxr-x--- 4 sunchengquan sunchengquan 4.0K 11月 12 2016 rmDuplicate
-rw-r----- 1 sunchengquan sunchengquan 103K 11月 12 2016 rmDuplicate.zip
-rw-r----- 1 sunchengquan sunchengquan 3.1K 5月 18 2017 test.bed
[sunchengquan 17:46:21 ~/test/sun]
$ tree rmDuplicate
rmDuplicate
├── picard
│ ├── paired
│ │ ├── readme.txt
│ │ ├── tmp.header
│ │ ├── tmp.MarkDuplicates.log
│ │ ├── tmp.metrics
│ │ ├── tmp.rmdup.bai
│ │ ├── tmp.rmdup.bam
│ │ ├── tmp.sam
│ │ └── tmp.sorted.bam
│ └── single
│ ├── readme.txt
│ ├── tmp.header
│ ├── tmp.MarkDuplicates.log
│ ├── tmp.metrics
│ ├── tmp.rmdup.bai
│ ├── tmp.rmdup.bam
│ ├── tmp.sam
│ └── tmp.sorted.bam
└── samtools
├── paired
│ ├── readme.txt
│ ├── tmp.header
│ ├── tmp.rmdup.bam
│ ├── tmp.rmdup.vcf.gz
│ ├── tmp.sam
│ ├── tmp.sorted.bam
│ └── tmp.sorted.vcf.gz
└── single
├── readme.txt
├── tmp.header
├── tmp.rmdup.bam
├── tmp.rmdup.vcf.gz
├── tmp.sam
├── tmp.sorted.bam
└── tmp.sorted.vcf.gz
6 directories, 30 files
十、打开第九题解压的文件,进入 rmDuplicate/samtools/single 文件夹里面,查看后缀为 .sam 的文件,搞清楚 生物信息学里面的SAM/BAM 定义是什么。
序列比对的存储格式,由一些比对软件产生,如bwa,bowtie2
sam stands for Sequence Alignment Mapping
sam‘ 序列比对映射’的首字母缩写
sam分为两部分,注释信息(header section)和比对结果部分(alignment section)
比对结果部分(alignment section),每一行表示一个片段(segment)的比对信息,包括11个必须的字段(mandatory fields)和一个可选的字段,字段之间用tag分割。必须的字段有11个,顺序固定,不可用时,根据字段定义,可以为’0‘或者’*’
- 1 QNAME,序列的名字(Read的名字)
- 2 FLAG, 概括出一个合适的标记,各个数字分别代表
- 3 RNAME,参考序列的名字(染色体)
- 4 POS,在参考序列上的位置(染色体上的位置)
- 5 MAPQ, mapping qulity 越高则位点越独特
- 6 CIGAR,代表比对结果的CIGAR字符串
- 7 RNEXT, mate 序列所在参考序列的名称; 下一个片段比对上的参考序列的编号,没有另外的片段,这里是’*‘,同一个片段,用’=‘;
- 8 PNEXT, mate 序列在参考序列上的位置;下一个片段比对上的位置,如果不可用,此处为0;
- 9 TLEN,估计出的片段的长度,当mate 序列位于本序列上游时该值为负值。Template的长度,最左边得为正,最右边的为负,中间的不用定义正负,不分区段(single-segment)的比对上,或者不可用时,此处为0
- 10 SEQ,read的序列;序列片段的序列信息,如果不存储此类信息,此处为’*‘,注意CIGAR中M/I/S/=/X对应数字的和要等于序列长度;
11 QUAL,ASCII码格式的序列质量;序列的质量信息,格式同FASTQ一样。
可选的字段(field) - 12 NM:i 经过编辑的序列
-13 MD:Z 代表序列和参考序列错配的字符串 - 14 AS:i 匹配的得分
bam是其压缩格式,samtools可以转化
samtools view -bS tmp.sam >tmp.bam
十一、安装 samtools 软件
curl -OL https://sourceforge.net/projects/samtools/files/samtools/1.6/samtools-1.6.tar.bz2
tar jxvf samtools-1.6.tar.bz2
cd samtools-1.6
./configure
make
mv samtools-1.6 samtools
ln -s ~/local/app/samtools/samtools ~/bin/samtools
查看samtools的手册
man ~/local/app/samtools/samtools.1
十二、打开 后缀为BAM 的文件,找到产生该文件的命令。 提示一下命令是:
/home/jianmingzeng/biosoft/bowtie/bowtie2-2.2.9/bowtie2-align-s --wrapper basic-0 -p 20 -x /home/jianmingzeng/reference/index/bowtie/hg38 -S /home/jianmingzeng/data/public/allMouse/alignment/WT_rep2_Input.sam -U /tmp/41440.unp
[sunchengquan 21:15:22 ~/test/sun/rmDuplicate/samtools/single]
$ samtools view -h tmp.sorted.bam |grep '^@PG'|awk 'BEGIN{FS="\t"}{print $5}'|cut -d: -f2
"/home/jianmingzeng/biosoft/bowtie/bowtie2-2.2.9/bowtie2-align-s --wrapper basic-0 -p 20 -x /home/jianmingzeng/reference/index/bowtie/hg38 -S /home/jianmingzeng/data/public/allMouse/alignment/WT_rep2_Input.sam -U /tmp/41440.unp"
十三题、根据上面的命令,找到我使用的参考基因组 /home/jianmingzeng/reference/index/bowtie/hg38 具体有多少条染色体
[sunchengquan 22:02:03 ~/test/sun/rmDuplicate/samtools/single]
$ samtools view -h tmp.sorted.bam |egrep '^@S.*?(chr[XYM]\s+.*|chr[1-9]?[0-9]\s+).*'|less
[sunchengquan 22:03:49 ~/test/sun/rmDuplicate/samtools/single]
$ samtools view -h tmp.sorted.bam |egrep '^@S.*?(chr[XYM]\s+.*|chr[1-9]?[0-9]\s+).*'|wc -l
25
十四题、上面的后缀为BAM 的文件的第二列,只有 0 和 16 两个数字,用 cut/sort/uniq等命令统计它们的个数
[sunchengquan 22:10:11 ~/test/sun/rmDuplicate/samtools/single]
$ samtools view tmp.rmdup.bam |cut -f2|sort |uniq -c
16 0
12 16
十五题、重新打开 rmDuplicate/samtools/paired 文件夹下面的后缀为BAM 的文件,再次查看第二列,并且统计
[sunchengquan 22:14:04 ~/test/sun/rmDuplicate/samtools/paired]
$ samtools view tmp.rmdup.bam |cut -f2 |sort |uniq -c|sort -t' ' -nrk1,1
8 99
7 147
2 97
2 83
2 163
1 433
1 387
1 371
1 353
1 323
十六题、下载 http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip 文件,并且解压,查看里面的文件夹结构, 这个文件有2.3M,注意留心下载时间及下载速度。
[sunchengquan 22:15:31 ~/test/sun]
$ wget -c http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip
[sunchengquan 22:21:36 ~/test/sun]
$ unzip sickle-results.zip
[sunchengquan 22:22:36 ~/test/sun]
$ tree sickle-results
sickle-results
├── command.txt
├── single_tmp_fastqc.html
├── single_tmp_fastqc.zip
├── test1_fastqc.html
├── test1_fastqc.zip
├── test2_fastqc.html
├── test2_fastqc.zip
├── trimmed_output_file1_fastqc.html
├── trimmed_output_file1_fastqc.zip
├── trimmed_output_file2_fastqc.html
└── trimmed_output_file2_fastqc.zip
0 directories, 11 files
十七题、解压 sickle-results/single_tmp_fastqc.zip 文件,并且进入解压后的文件夹,找到 fastqc_data.txt 文件,并且搜索该文本文件以 >>开头的有多少行?
[sunchengquan 22:23:20 ~/test/sun/sickle-results]
$ unzip single_tmp_fastqc.zip
[sunchengquan 22:23:42 ~/test/sun/sickle-results]
$ cd single_tmp_fastqc
[sunchengquan 22:24:06 ~/test/sun/sickle-results/single_tmp_fastqc]
$ ls
fastqc_data.txt fastqc.fo fastqc_report.html Icons Images summary.txt
[sunchengquan 22:25:02 ~/test/sun/sickle-results/single_tmp_fastqc]
$ grep '^>>' fastqc_data.txt |wc -l
24
[sunchengquan 08:02:54 ~/test/sun/sickle-results/single_tmp_fastqc]
$ cat fastqc_data.txt |awk '/^>>/{print}'|wc -l
24
十八题、下载 http://www.biotrainee.com/jmzeng/tmp/hg38.tss 文件,去NCBI找到TP53/BRCA1等自己感兴趣的基因对应的 refseq数据库 ID,然后找到它们的hg38.tss 文件的哪一行。
[sunchengquan 22:27:31 ~/test/sun]
$ wget -c http://www.biotrainee.com/jmzeng/tmp/hg38.tss
[sunchengquan 22:43:17 ~/test/sun]
$ grep 'NM_000546' hg38.tss
NM_000546 chr17 7685550 7689550 1
[sunchengquan 22:43:27 ~/test/sun]
$ grep 'NM_001126113' hg38.tss
NM_001126113 chr17 7685550 7689550 1
十九题、解析hg38.tss 文件,统计每条染色体的基因个数。
[sunchengquan 22:50:12 ~/test/sun]
$ cat hg38.tss |cut -f2|sort|uniq -c
6050 chr1
2824 chr10
………
[sunchengquan 22:59:12 ~/test/sun]
$ grep -oE 'chr[0-9]{1,2}|chr[a-zA-Z]{1,2}' hg38.tss |sort |uniq -c
6157 chr1
2838 chr10
3577 chr11
3014 chr12
1133 chr13
1982 chr14
2377 chr15
2696 chr16
3794 chr17
883 chr18
5880 chr19
4090 chr2
1692 chr20
895 chr21
1410 chr22
3395 chr3
2277 chr4
2821 chr5
5782 chr6
2785 chr7
2221 chr8
2310 chr9
2 chrM
32 chrUn
2561 chrX
414 chrY
二十题、解析hg38.tss 文件,统计NM和NR开头的序列,了解NM和NR开头的含义。
[sunchengquan 23:01:06 ~/test/sun]
$ grep '^NR' hg38.tss |wc -l
15954
[sunchengquan 23:01:13 ~/test/sun]
$ grep '^NM' hg38.tss |wc -l
51064
$ grep -oE '^(NM|NR)' hg38.tss |sort|uniq -c
51064 NM
15954 NR
[sunchengquan 23:06:29 ~/test/sun]
$ grep -E '^(NM|NR)' hg38.tss |wc -l
67018
[sunchengquan 23:08:30 ~/test/sun]
$ wc -l hg38.tss
67018 hg38.tss
NM:转录组产物的序列mRNA
NR:非编码的转录组序列ncRNA