一到十二题由于比较简单,所以我直接给出了命令和输出。十二题后详述了解决方法和做题思路。这次做题后对grep、cut 、awk、wc等命令印象更加深刻,学到很多知识。但是学生初来乍到,有错误在所难免,还请老师们纠正。
以下是我完成得作业:
一、在任意文件夹下面创建形如 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 ,内容也要一样!
7.png
七,再次删除掉前面几个步骤建立的文件夹及文件
八、下载 http://www.biotrainee.com/jmzeng/igv/test.bed 文件,后在里面选择含有 H3K4me3 的那一行是第几行,该文件总共有几行。
九、下载 http://www.biotrainee.com/jmzeng/rmDuplicate.zip 文件,并且解压,查看里面的文件夹结构
解压:unzip
sam ,bam文件详解点这里,这个作者把这个讲得很详细:
https://pzweuj.github.io/2019/02/15/BAM.html
十、打开第九题解压的文件,进入 rmDuplicate/samtools/single 文件夹里面,查看后缀为 .sam 的文件,搞清楚 生物信息学里面的SAM/BAM 定义是什么。
十一、安装 samtools 软件:
安装软件conda install samtools
十二、打开 后缀为BAM 的文件,找到产生该文件的命令。
十三题、根据上面的命令,找到我使用的参考基因组具体有多少条染色体
grep 'chr[1-9,X,Y]' /teach/database/genome/hg38.fa : 查看hg38文件中染色体列
cut -f1|sort|uniq|wc -l:提取染色体列并进行计数
grep 'chr[1-9,X,Y]' /teach/database/genome/hg38.fa |cut -f1|sort|uniq|wc -l
命令效果如下,这样就一次性达到了老师要求。我之前查看了其他同学的作业,需要输出好几个文本,挨个去查看和计数,因为我对于取文件的名字实在是头痛,就直接一次性输出统计数了。管道符大家一定要用起来,会起到事半功倍的效果。
wc命令的具体使用方法点这里:https://www.cnblogs.com/peida/archive/2012/12/18/2822758.html
十四题、上面的后缀为BAM 的文件的第二列,只有 0 和 16 两个数字,用 cut/sort/uniq等命令统计它们的个数。
开始想直接用cut命令切出文本,结果一串乱码出来。bam文件是sam 文件的二进制格式,查看需要用到samtool view命令。用这个命令将需要的统计的文本打开,然后再用cut命令切出文本第二行,grep 命令查看指定字段,wc 统计行数,大功告成!
1,查看0的个数
samtools view tmp.rmdup.bam|cut -f 2 |grep '0'|wc -l
2,查看16的个数
samtools view tmp.rmdup.bam|cut -f 2 |grep '16'|wc -l
十五题、重新打开 rmDuplicate/samtools/paired 文件夹下面的后缀为BAM 的文件,再次查看第二列,并且统计
samtools view tmp.rmdup.bam|cut -f 2|sort|uniq -c
uniq 命令参数用法看这里:
https://www.runoob.com/linux/linux-comm-uniq.html
注:当重复的行并不相邻时,uniq 命令是不起作用的,即若文件内容为以下时,uniq 命令不起作用,这就是为什么先用sort的原因。
十六题、下载 http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip 文件,并且解压,查看里面的文件夹结构, 这个文件有2.3M,注意留心下载时间及下载速度。
1,下载
2,解压
3,tree 命令查看文件夹整体结构
十七题、解压 sickle-results/single_tmp_fastqc.zip 文件,并且进入解压后的文件夹,找到 fastqc_data.txt 文件,并且搜索该文本文件以 >>开头的有多少行?
十八题、下载 http://www.biotrainee.com/jmzeng/tmp/hg38.tss 文件,去NCBI找到TP53/BRCA1等自己感兴趣的基因对应的 refseq数据库 ID,然后找到它们的hg38.tss 文件的哪一行。
1,登录ncbi,选择gene,找到下图部分,可以看到每个剪切变体前的以NM开头的为id号。
2,可以用grep,或者cut 命令查看:
十九题、解析hg38.tss 文件,统计每条染色体的基因个数。
tss 文件内容如下:
处理思路和第十三题有些相似:可以由上观察到,每个染色体号名字后面有多个基因,有的不同基因对应着相同的染色体号。所以可以先把此文件第2列提取出来,然后用sort和uniq命令对其进行归类和统计。对文本中指定内容进行提取可以用cut,或者awk,故本题解决方法可以有两种:
1,使用cut
cut -f2 hg38.tss|sort|uniq -c|head -n10
输出结果如下:
2,使用awk
awk '{print $2}' hg38.tss |sort |uniq -c
二十题、解析hg38.tss 文件,统计NM和NR开头,了解NM和NR开头的含义。
1,按指定字段查看hg38.tss 文件
2,对NM和NR进行统计
用指定字符"_",对文本进行分割,提取文本NM和NR开头,再用|sort|uniq -c进行统计
cut -d "_" -f1 hg38.tss |head -n10
效果如下:
NM:mRNA mixed,转录组产物序列;成熟mRNA转录本序列
NR:RNA mixed,非编码的转录子序列,包括结构RNAs,假基因转子等