一、在任意文件夹下面创建形如 1/2/3/4/5/6/7/8/9 格式的文件夹系列。

mkdir –p 1/2/3/4/5/6/7/8/9

pwd

二、在创建好的文件夹下面，比如我的是 /tmp/test/1/2/3/4/5/6/7/8/9 ，里面创建文本文件 me.txt。

touch me.txt

三、在文本文件 me.txt 里面输入内容:

Go to: http://www.biotrainee.com/

I love bioinfomatics.

And you ?

cat > me.txt

Go to: http://www.biotrainee.com/

I love bioinfomatics.

And you ?

Ctrl+d

四、删除上面创建的文件夹 1/2/3/4/5/6/7/8/9 及文本文件 me.txt

rm –rf 1/

rm –rf me.txt

五、在任意文件夹下面创建 folder1~5这5个文件夹，然后每个文件夹下面继续创建 folder1~5这5个文件夹

mkdir –p folder_{1..5}/folder_{1..5}

六、在第五题创建的每一个文件夹下面都 创建第二题文本文件 me.txt ，内容也要一样。

touch folder_{1..5}/folder_{1..5}/me.txt

vim me.txt (输入内容)

echo folder_{1..5}/ folder_{1..5}|xargs –n 1

echo folder_{1..5}/ folder_{1..5}|xargs –n 1 cp me.txt –v

ls * / * /

七、再次删除掉前面几个步骤建立的文件夹及文件

rm –rf folder_*

rm –rf me.txt

wget –c http://www.biotrainee.com/jmzeng/igv/test.bed

ls

pwd

grep -n -o --color H3K4me3 /root/test/test.bed

(-n 标记行数，-o 只显示匹配上的，--color匹配文字出现颜色)

cat test.bed |wc -l

(wc显示文件的行数、单词数、字节数)

wget -c http://www.biotrainee.com/jmzeng/rmDuplicate.zip

ls

yum install -y unzip zip；

unzip rmDuplicate.zip

yum install -y tree

ls

tree rmDuplicate

十、打开第九题解压的文件，进入rmDuplicate/samtools/single 文件夹里面，查看后缀为 .sam 的文件，搞清楚 生物信息学里面的SAM/BAM 定义是什么。

cd rmDuplicate/samtools/single

ls

cat tmp.sam | head

SAM/BAM 定义是什么

Sam:头部注释信息和比对结果，每行包含12列，以tab分隔

Bam：SAM的二进制文件，内存小

1)注释信息可有可无，都是以@开头

@HD，说明符合标准的版本、对比序列的排列顺序；

@SQ，参考序列说明；

@RG，比对上的序列(read)说明；

@PG，使用的程序说明；

@CO，任意的说明信息。

2)详细比对：11个tab隔开的字段

第一列：序列名称

第二列：序列标记信息

1 序列是一对序列中的一个

2 比对结果是一个pair-end比对的末端

4 没有找到位点

8 这个序列是pair中的一个但是没有找到位点

16 在这个比对上的位点，序列与参考序列反向互补

32 这个序列在pair-end中的的mate序列与参考序列反响互补

64 序列是 mate 1

128 序列是 mate 2

如果以上情况都不符合则默认为0

第三列：比对上的参考基因染色体号。没比对上则显示*

第四列：比对上的在参考基因组上的起始位置，没有比对上则为0

第五列：比对质量分数。数字越大，特异性越高。但值为255是指比对质量不可用

第六列：比对结果的CIGAR 的字符串

(即比对的详细情况， 记录插入，删除，错配，后剪切拼接的接头 。以参考序列为基础，使用数字+字母表示比对结果。)

单端测序没有第7-9

第七列：RNEXT ，双端测序中下一个reads比对的参考系列的名称。“＊”是完全没有比对上，“＝”代表完全比对

第八列：PNEXT 如果是双端测序，是指另一端匹配到参考基因组的位置，如果设置为0，那么该列不可用

第九列： TLEN Template的长度，最左边得为正，最右边的为负，中间的不用定义正负，不分区段(single-segment)的比对上，或者不可用时，此处为0；

第十列：SEQ序列片段的序列信息，如果不存储此类信息，此处为’*‘，注意CIGAR中M/I/S/=/X对应数字的和要等于序列长度；

第十一列：QUAL序列的质量信息，格式同FASTQ一样。ASCII码

第十二列：可选字段 格式：TAG:TYPE:VALUE TAG是两个大写字母，每个TAG代表一类信息。TYPE代表TAG对应值的类型(字符串、数组、字节等)

十一、安装 samtools 软件

(提前预装好conda)

conda install samtools

输入`y

十二、打开 后缀为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

cd rmDuplicate/samtools/single

ls

which samtools(查找路径)

samtools(路径) view -H tmp.rmdup.bam

产生bam文件的命令一般都在头部注释信息中，但都会放在最后用tail

tail -n 3 tmp.header

十三题、根据上面的命令，找到我使用的参考基因组 /home/jianmingzeng/reference/index/bowtie/hg38 具体有多少条染色体。

samtools view -H tmp.rmdup.bam |grep -o -E "SN:chr[0-9]+|SN:chr[a-zA-Z]+"|sort |uniq -c |wc -l

十四题、上面的后缀为BAM 的文件的第二列，只有 0 和 16 两个数字，用 cut/sort/uniq等命令统计它们的个数。

samtools(路径) view tmp.rmdup.bam | cut -f 2 | sort -n | uniq -dc

十五题、重新打开 rmDuplicate/samtools/paired 文件夹下面的后缀为BAM 的文件，再次查看第二列，并且统计

cd rmDuplicate/samtools/paired

ls

samtools(路径) view tmp.rmdup.bam | cut -f 2 | sort -n | uniq -dc

十六题、下载 http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip 文件，并且解压，查看里面的文件夹结构， 这个文件有2.3M，注意留心下载时间及下载速度。

十七题、解压 sickle-results/single_tmp_fastqc.zip 文件，并且进入解压后的文件夹，找到 fastqc_data.txt 文件，并且搜索该文本文件以 >>开头的有多少行？

unzip sickle-results/single_tmp_fastqc.zip

cd single_tmp_fastqc

cat fastqc_data.txt | grep '>>' | wc -l

十八题、下载 http://www.biotrainee.com/jmzeng/tmp/hg38.tss 文件，去NCBI找到TP53/BRCA1等自己感兴趣的基因对应的 refseq数据库 ID，然后找到它们的hg38.tss 文件的哪一行。(https://www.ncbi.nlm.nih.gov/gene/7157)

wget -c http://www.biotrainee.com/jmzeng/tmp/hg38.tss

grep NM_000546 hg38.tss

关于它的ID：NM开头的表示标准序列，XM表示预测的蛋白编码序列，NR表示非编码蛋白的mRNA序列，AF开头的表示克隆序列，BC开头的表示模板序列

另外，你可能见过gi|4557284|ref|NM_000646.1|[4557284]这种格式

gi就是代表genebank identifier；ref就是对应的refseq中的ID啦

十九题、解析hg38.tss 文件，统计每条染色体的基因个数。(2种方法)

less hg38.tss

cut -f2 hg38.tss |cut -d'_' -f1 |sort |uniq -c |sort -rn

cat hg38.tss |cut -f 2|grep -o -E "chr[0-9]{1,2}|chr[a-zA-Z]{1,2}"|sort |uniq -dc

二十题、解析hg38.tss 文件，统计NM和NR开头的序列，了解NM和NR开头的含义。(3种方法)

cat hg38.tss |grep 'NM' | wc -l

cat hg38.tss |grep 'NR' | wc -l

grep -o -E "NM|NR" hg38.tss |sort |uniq -dc

cat hg38.tss |awk '{print $1}'|cut -c1-2|sort|uniq -c

(只带有前2个字母进行统计)

NM:开头表示标准序列。可以转录成蛋白质的基因

NR:开头表示非编码蛋白的mRNA序列

数据练习

(首先下载好数据)

看文件大小

ll -h *.gff

看文本内容

less -SN *.gff

去除多余的#注释行与 空行，并打印行号

grep -v "#" *.gff | grep -v "^$" | wc -l

截取文件1-5列，将第2列去除，输出默认的前10行到test.txt中

cut -f 1,3,4,5 *.gff > test.txt

cat test.txt | head

(然后对截取的test.txt进行处理)

根据第2列的feature进行排序

sort -k 2,2 test.txt | head -n5

先根据第1列Chr数字大小降序排序(第一个字段的第四个字符)，再根据第2列排序

sort -k 1.4,1.4nr -k 2,2 *.txt | head -n5

统计整体的feature(第3列)

cut -f 3 *.gff |sort|uniq -c >feature.txt

哪个feature最多

sort -k 1r feature.txt

练习：

GENCODE下载人类基因组GRCh38注释gff3，然后统计人类基因组feature(第1列 参照序列，2来源，3类型，4起点，5终点，6得分，7链(正负链)，8步进，9属性)

1.下载

wget -c ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz

2.解压

gunzip gencode.v28.annotation.gff3.gz

3.看文本内容

less -SN *.gff3

4.去除多余的#注释行与 空行，统计整体的feature(第3列)

grep -v "#" *.gff3 | grep -v "^$" | cut -f 3 |sort|uniq -c >feature.txt

5.查看文件(前10行)

cat feature.txt | head

哪个feature最少(从小到大)

sort -k 1n feature.txt

(从大到小)

sort -k1,1nr feature.txt